#!/bin/perl
#
# 1 pixel 当たりの発生BG electronを見積もる
# そうすると
#
# 背景放射密度    I_nu_Mjy_str   MJy/str
# 波長 lambda (m) lambda_mu      um            
# 口径 (m)        D              m
# 比波長分解能    R  
# 効率 eta (e/ph) eta
# 積分時間        T              s
#
#
# pixel scale 回折限界をナイキストサンプル
# 回折限界 lamda/D (rad)
# x_pix=lamda/D/2  (rad)
#
# lambda=lambda_mu * 1e-6
#
# I_lambda_BG = I_nu_BG*d_nu/d_lamda = I_nu_BG*c**2/lamda
# Omega_pixel = (x_pix)**2 = (lamda/2D)**2
# A_tel       = pi*(D/2)**2
# E_photon    = hc/lamda
# d_lamda     = R/lamda

# singal lamda background
#   *=  I_lambda_BG          W m-2 m-1 str-1
#   *=  Omega_pixel          W m-2 m-1
#   *=  d_lamda              W m-2
#   *=  A_telescope          W
#   /=  E_photon             photon/s
#   *=  eta                  e/s
use strict;

# flag
my $debug=1;
my $latex=1;

# constants
my $pi=4*atan2(1,1);         # pai
my $c=3.00e8;                # speed of light m/s
my $h=6.63e-34;              # Planck constant Js

# printf(STDERR "# pi %e\n",$pi) if($debug);

#   um    MJy/str
#   12    10.7
#   25    18
#  100    2.5
#printf("%20s %20s %20s %20s\n","wavelength(um)","sky brightness MJy/str","sky signal electron/s/pixel","pixel scale (arcsec)");
if($latex){
    printf("%15s &%15s &%15s &%15s\\\\\n","wavelength","sky brightness","sky signal","pixel scale");
    printf("%15s &%15s &%15s &%15s\\\\\n","(um)","MJy/str","e/s/pixel","arcsec");
}else{
    printf("%15s %15s %15s %15s\n","wavelength","sky brightness","sky signal","pixel scale");
    printf("%15s %15s %15s %15s\n","um","MJy/str","e/s/pixel","arcsec");
}

&fir_bg_eletrons(1.25,0.16);
&fir_bg_eletrons(2.2,0.11);
&fir_bg_eletrons(3.5,0.074);
&fir_bg_eletrons(4.9,0.31);
&fir_bg_eletrons(12,10.7);
&fir_bg_eletrons(25,18);
&fir_bg_eletrons(60,5.2);
&fir_bg_eletrons(100,2.5);
&fir_bg_eletrons(140,2.7);
&fir_bg_eletrons(240,1.8);
 
sub fir_bg_eletrons{
# parameters
#my $I_nu_MJy    =2.5;    # MJy/str
#my $lambda_mu    =100;    # micron
my $lambda_mu    =$_[0];
my $I_nu_MJy     =$_[1];

my $D           =3.0;   # m
#my $D           =1.0;   # m
my $R           =10;    # 
my $eta         =0.15;   # electron/s  
my $T           =1;     # s
#
printf(STDERR "# parameters\n") if($debug);
printf(STDERR "# background          %g MJy/str\n",$I_nu_MJy) if($debug);
printf(STDERR "# wavelength          %g micron\n",$lambda_mu) if($debug);
printf(STDERR "# telescope diameter  %g m\n",$D) if($debug);
printf(STDERR "# spectral resolution %g\n",$R) if($debug);
printf(STDERR "# efficiency          %g e/photon\n",$eta) if($debug);
#
my $I_nu = $I_nu_MJy*1e6*1e-26;  # MJy/str => W/m2/Hz/str
my $lambda= $lambda_mu*1e-6     ;  # micron  => m

# 
my $I_lambda = $I_nu*$c/$lambda**2;  # W/m2/Hz => W/m2/m
my $x_pix=$lambda/$D/2;              # pixel scale, nyquist of difraction limit
my $Omega_pixel = ($x_pix)**2;      # 立体角 radian => steratian
my $A_telescope = $pi*($D/2)**2;    # 有効面積
my $E_photon    = $h*$c/$lambda;      # photon energy (J)
my $d_lambda     = $lambda/$R;        # 帯域

printf(STDERR "# I_lamda %e W/m2/str\n",$I_lambda) if($debug);

my $x_pix_arcsec=$x_pix/4.85e-6;    # radian  =>  arcsec
printf(STDERR "# pixel scape %g arcsec\n",$x_pix_arcsec) if($debug);
my $nu=$c/$lambda;                  # freq

my $singal_lambda_bg=1;
$singal_lambda_bg*=$I_lambda;
printf(STDERR "# I_lambda %g W/m2 m str\n",$singal_lambda_bg) if($debug);
$singal_lambda_bg*=$Omega_pixel;
printf(STDERR "# I_lambda Omega_pix %g W/m2 m pix\n",$singal_lambda_bg) if($debug);
$singal_lambda_bg*=$d_lambda;
printf(STDERR "# I_lambda Omega_pix d_lambda %g W/m2 pix\n",$singal_lambda_bg) if($debug);
$singal_lambda_bg*=$A_telescope;
printf(STDERR "# I_lambda Omega_pix d_lambda A_tel %g W/pix\n"
       ,$singal_lambda_bg) if($debug);
$singal_lambda_bg/=$E_photon;
printf(STDERR "# I_lambda Omega_pix d_lambda A_tel (E_photon)-1  %g photon/s pix\n"
       ,$singal_lambda_bg) if($debug);
$singal_lambda_bg*=$eta;
printf(STDERR "# I_lambda Omega_pix d_lambda A_tel (E_photon)-1 eta %g electron/s pix\n"
       ,$singal_lambda_bg) if($debug);
printf(STDERR "# I_lambda Omega_pix d_lambda A_tel (E_photon)-1 eta %g electron/s pix\n"
       ,$I_lambda*($lambda**4)/16*$pi/($h*$c)*$eta/$R) if($debug);

if($latex){
    printf("%15g &%15g &%15.1e &%15.2g\\\\\n",$lambda_mu,$I_nu_MJy,$singal_lambda_bg,$x_pix_arcsec);
}else{
    printf("%15g %15g %15.1e %15.2g\n",$lambda_mu,$I_nu_MJy,$singal_lambda_bg,$x_pix_arcsec);
}

}