#! /usr/bin/perl -T

## Calculator for Limiting Magnitude
## Use at your own risk. There is no guarantee, nor supprt for this software.
## 2008/9/3 Kentaro Motohara (kmotohara@ioa.s.u-tokyo.ac.jp)

#物理定数
$C=3e8;
$H=6.626e-34;
$K=1.38e-23;
$PI=3.1415926;

$narg = 4;
if ( $#ARGV != ($narg - 1) ) {
print("Usage limmag (D : m) (exposure : sec) (S/N)\n");
exit(1);
}
$diameter=$ARGV[0];
$exposure=$ARGV[1];
$sn=$ARGV[2];
$paramfile=$ARGV[3];

open(INSTFILE,"< instrument.dat");
while (<INSTFILE>){
    @line= split;
    if($line[0] eq "Dark"){$dark_opt=$line[1];$dark_nir=$line[2];$dark_mir=$line[3];}
    if($line[0] eq "Well"){$well_opt=$line[1];$well_nir=$line[2];$well_mir=$line[3];}
    if($line[0] eq "RN"){$readnoise_opt=$line[1];$readnoise_nir=$line[2];$readnoise_mir=$line[3];}
    if($line[0] eq "ExpMin"){$sexpmin_opt=$line[1];$sexpmin_nir=$line[2];$sexpmin_mir=$line[3];}
}
close(INSTFILE);

# Read Filter Parameters from $paramfile
open(PARAMFILE,"< $paramfile");
while (<PARAMFILE>){
    @line= split;
    if($line[0] eq "Name"){shift(@line);@filter=@line}
    elsif($line[0] eq "Lambda"){shift(@line);@lambda=@line}
    elsif($line[0] eq "DLambda"){shift(@line);@dlambda=@line}
    elsif($line[0] eq "ZMFlux"){shift(@line);@zmflux=@line}
    elsif($line[0] eq "BGMag"){shift(@line);@bgmag=@line}
    elsif($line[0] eq "Nseeing"){shift(@line);@seeing_natural=@line}
    elsif($line[0] eq "AOcore"){shift(@line);@imagesize=@line}
    elsif($line[0] eq "PixScale"){shift(@line);@pixscale=@line}
    elsif($line[0] eq "Aperture"){shift(@line);@diam=@line}
    elsif($line[0] eq "SR"){shift(@line);@strehl=@line}
    elsif($line[0] eq "Eatm"){shift(@line);@thrputatm=@line}
    elsif($line[0] eq "Etel"){shift(@line);@thrputtel=@line}
    elsif($line[0] eq "Einst"){shift(@line);@thrputinst=@line}
    elsif($line[0] eq "Eao"){shift(@line);@thrputao=@line}
}
close(PARAMFILE);

# If throughput == 0, set it to 1e-50
# 大気の透過率が20%以下の領域は観測できないとする。($thrputmisc[$i]=1e-50に)。
for( $i=0; $i<$#filter+1 ; $i++ ){
    $thrputmisc[$i]=1.0; # まずはthrputmiscすべてが1だとする。
    if($thrputatm[$i] == 0) { $thrputatm[$i] = 1e-50; }
    if( $thrputatm[$i] < 0.2 ){	$thrputmisc[$i]=1e-80; } 
}



# Calculate Aperture Size
for( $i=0; $i<$#filter+1 ; $i++ ){
    $aperture[$i]=int($diam[$i]*$diam[$i]/4*$PI)+1;
}

# Calculate Nu
for( $i=0; $i<$#filter+1 ; $i++ ){
    $nu[$i]=$C/$lambda[$i];
}

# Calculate size of diffraction limit at each band
# 2.44*lam/D
for( $i=0; $i<$#filter+1 ; $i++){
    $seeing_diffrac[$i]=2.44*$lambda[$i]/$diameter/$PI*180*3600;
# If the seeing size is smaller than the diffraction limit, use diffraction limit size instead.
    if($seeing_diffrac[$i]>$seeing_natural[$i]){$seeing_natural[$i]=$seeing_diffrac[$i]}
# If $imagesize is "d", use diffraction limit size for it.
    if($imagesize[$i] eq "d"){$imagesize[$i]=$seeing_diffrac[$i]}
}

#Set throughput of Telescope+(AO)
for( $i=0; $i<$#filter+1 ; $i++ ){
    if ($strehl[$i]>0){$thrput[$i]=$thrputao[$i]*$thrputtel[$i]}
    else {$thrput[$i]=$thrputtel[$i]}
}

print("# D=$diameter(m)\n");
print("# T=$exposure(s)\n");
print("# S/N=$sn\n");
print("# Opt: Dark=$dark_opt (e-/s/pix)\tWell=$well_opt (e-/pix)\tReadnoise=$readnoise_opt (e- r.m.s/pix)\tMin_exposure=$sexpmin_opt (s)\n");
print("# Nir: Dark=$dark_nir (e-/s/pix)\tWell=$well_nir (e-/pix)\tReadnoise=$readnoise_nir (e- r.m.s/pix)\tMin_exposure=$sexpmin_nir (s)\n");
print("# Mir: Dark=$dark_mir (e-/s/pix)\tWell=$well_mir (e-/pix)\tReadnoise=$readnoise_mir (e- r.m.s/pix)\tMin_exposure=$sexpmin_mir (s)\n");
print("#\n");

for( $i=0; $i<$#filter+1 ; $i++ ){
    if($lambda[$i]<1e-6){$imagesize[$i]=$seeing_natural[$i];}
    else{$imagesize[$i]=$seeing_natural[$i];}
    $encirce[$i]=&encircled_e2($strehl[$i],$imagesize[$i],$seeing_natural[$i],$diam[$i]);
#	print("$seeing_diffrac[$i]\t$seeing_natural[$i]\t$imagesize[$i]\t$encirce[$i]\n")
}

#Set Readout Noise
for( $i=0; $i<$#filter+1 ; $i++ ){
    if($lambda[$i]<1e-6){$readnoise[$i]=$readnoise_opt;}
    elsif($lambda[$i]<5e-6){$readnoise[$i]=$readnoise_nir;}
    else{$readnoise[$i]=$readnoise_mir;}
}

#Set Dark Current
for( $i=0; $i<$#filter+1 ; $i++ ){
    if($lambda[$i]<1e-6){$dark[$i]=$dark_opt;}
    elsif($lambda[$i]<5e-6){$dark[$i]=$dark_nir;}
    else{$dark[$i]=$dark_mir;}
}

#背景放射による単位時間当りに検出される光子の数
for( $i=0; $i<$#filter+1 ; $i++ ){
#まずはフラックス
    $f=$zmflux[$i]*(exp(log(10.0)*(-0.4*$bgmag[$i])));
#検出されるPhoton数(/pix/sec)
    $ssky[$i]=$f*$diameter*$diameter/4*$PI*$dlambda[$i]*1e6*$thrput[$i]*$thrputinst[$i]/$H/$nu[$i]*$pixscale[$i]*$pixscale[$i];
}

#単一積分時間は
for( $i=0; $i<$#filter+1 ; $i++ ){
    if($ssky[$i]==0) { printf("%d %f\n",$i,$ssky[$i]); }
#まずはどのくらいでwellがいっぱいになるか
#もしもsexpminだとwellが一杯になってしまうようだと、その波長の観測はできないことにする($thrputmisc[$i]=1e-60に)。
    if($lambda[$i]<1e-6){
# 可視では5分
	$sexposure[$i]=300
    }
    elsif($lambda[$i]<2.5e-6){
# キンセキガイでは2分
	$sexposure[$i]=$well_nir/$ssky[$i];
	if($sexposure[$i]<$sexpmin_nir){$thrputmisc[$i]=1e-60;}
    }
    elsif($lambda[$i]<5e-6){
	$sexposure[$i]=$well_nir/$ssky[$i];
	if($sexposure[$i]<$sexpmin_nir){$thrputmisc[$i]=1e-60;}
    }
    else{
	$sexposure[$i]=$well_mir/$ssky[$i];
	if($sexposure[$i]<$sexpmin_mir){$thrputmisc[$i]=1e-60;}
    }

#sexposureが3600秒を越えるのであれば強制的に3600秒にしてしまう。
    if($sexposure[$i] > 3600) {
	$sexposure[$i]=3600;
    }

#つぎに、どのくらいでbackground limitedになるか
# 4 Readnoise = BGnoise が分岐点で、それに到達する積分時間をminexposureとする

    $minexposure[$i]=16*$readnoise[$i]*$readnoise[$i]/$ssky[$i];

#minexposureが20秒以下であれば強制的に20秒にしてしまう。
    if($minexposure[$i] < 20) {
	$minexposure[$i]=20;
    }
#    print("$minexposure[$i] ");

##minexposure と sexposure を比較して、小さい方を採用する。
    if( $minexposure[$i] < $sexposure[$i] ){
	$sexposure[$i]=$minexposure[$i];
    }

#総積分時間でwellがあふれなければ…
    if($sexposure[$i] >= $exposure) {
	$sexposure[$i]=$exposure;
	$n[$i]=1;
    }
#総積分時間でwellがあふれるなら…
    else{
	$n[$i]=int($exposure/$sexposure[$i]);
	if($n[$i] != $exposure/$sexposure[$i]) {$n[$i]+=1;}
	$sexposure[$i]=$exposure/$n[$i];
    }
#    printf("%d %.2e %d\t %.3f %.3f\n",$i,$pixscale[$i],$n[$i],$minexposure[$i],$sexposure[$i]);
}

printf("# Lam      dLam    F0_nu(Jy)\tExp      x  Number\tF_lam    F_nu(Jy) Mag\tPixscale Bin\tS_sky    S_obj\t\tImagesize\tStrehl\tEatm\tEinst\tEtel\tEmisc\n");
for( $i=0; $i<$#filter+1 ; $i++ ){
    $noise[$i]=$aperture[$i]*($ssky[$i]+$dark[$i]+$readnoise[$i]*$readnoise[$i]/$sexposure[$i]);
#    printf("%e %e %d\n",$noise[$i],$ssky[$i]*$aperture[$i],$aperture[$i]);
   $sobj[$i]=(($sn**2)+$sn*sqrt(($sn**2)+4*$exposure*$noise[$i]))/2/$exposure;
#Encircled Energy の分を戻す。
    $sobj[$i]/=$encirce[$i];
#    $sobj[$i]=$sn/sqrt($exposure)*sqrt($noise[$i]);

#  noiseother=ssky+sbb+pixscalebin*pixscalebin*dark+pixscalebin*pixscalebin*readnoise*readnoise/tsample;
#  sobj=(sn*sn+sqrt(sn*sn*sn*sn+4*tint*sn*sn*noiseother))/(2*tint);
#    $flux[$i]=$sobj[$i]/($diameter**2)*4/$PI/$lambda[$i]/1e6*$lambda[$i]/$dlambda[$i]/$thrput[$i]*$H*$nu[$i];
    $flux[$i]=$sobj[$i]/$diameter/$diameter*4/$PI/1e6/$dlambda[$i]/$thrput[$i]/$thrputinst[$i]/$thrputatm[$i]/$thrputmisc[$i]*$H*$nu[$i];
    $fluxjy[$i]=$flux[$i]*1e6*$lambda[$i]*$lambda[$i]/$C*1e26;
    $zmfluxjy[$i]=$zmflux[$i]*1e6*$lambda[$i]*$lambda[$i]/$C*1e26;
    $mag[$i]=-2.5*log($flux[$i]/$zmflux[$i])/log(10.0);
    printf("%s %.3e %.2e %.2e\t%.2e x %7d\t%.2e %.2e %.2f\t%.2e %.1f\t%.2e %.2e\t%.2e\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\n",$filter[$i],$lambda[$i],$dlambda[$i],$zmfluxjy[$i],$sexposure[$i],$n[$i],$flux[$i],$fluxjy[$i],$mag[$i],$pixscale[$i],$aperture[$i],$ssky[$i],$sobj[$i],$imagesize[$i],$strehl[$i],$thrputatm[$i],$thrputinst[$i],$thrput[$i],$thrputmisc[$i]);
}




#Strehl Ratio(SR)によって、どの程度の光がもとめようとしているアパーチャーに入ってくるのかを計算。
# (完全なdiffraction limit)=(diffraction core)+(seeing halo)
# ∫exp(-r^2/r_0^2) dr^2= SR*∫exp(-r^2/r_0^2) dr^2 + A*∫exp(-r^2/r_1^2)dr^2
#   θ_diffrac=2.3548 * r_0
#   θ_seeing=2.3548 * r_1
#   A = シーイングハロのピーク
# θとSRが分かっていて、Aを求め、さらにあるアパーチャー内に入っているエネルギーの割合を算出。
# まずはAを求めるサブルーチン
sub halo_peak {
# $_[0]=strehl ratio
# $_[1]=r_0
# $_[2]=r_1
    local($lSR)=$_[0];
    local($lr0)=$_[1];
    local($lr1)=$_[2];
# ∫exp(-r^2/r_0^2) dr^2= 2π∫exp(-r^2/r_0^2) r dr=π * r_0^2 より
# A=(1-s)(r_0/r_1)^2
    return((1-$lSR)*$lr0*$lr0/$lr1/$lr1);
}

# 次にSR,A,θ_0, θ_1から、θの直径のencircled energy ratioをもとめる
sub encircled_e {
# $_[0]=strehl ratio
# $_[1]=A
# $_[2]=theta_0
# $_[3]=theta_1
# $_[3]=theta
    local($lSR)=$_[0];
    local($lA)=$_[1];
    local($lr0)=$_[2]/2.355;
    local($lr1)=$_[3]/2.355;
    local($lr)=$_[4]/2;
# Encircled Energy Ratio= (SR ∫_0^r exp(-r^2/r_0^2) dr^2 + A ∫_0^r exp(-r^2/r_1^2) dr^2 )/∫exp(-r^2/r_0^2) dr^2 
#                       = SR (1-exp(-r^2/r_0^2))+ A*r_1^2/r_0^2*(1-exp(-r^2/r_1^2))
    return($lSR*(1-exp(-$lr*$lr/$lr0/$lr0)) + $lA * $lr1*$lr1/$lr0/$lr0*(1-exp(-$lr*$lr/$lr1/$lr1)));
}

#上二つをまとめたもの
sub encircled_e2 {
# $_[0]=strehl ratio
# $_[2]=theta_0
# $_[3]=theta_1
# $_[4]=theta
    local($lSR)=$_[0];
    local($lr0)=$_[1]/2.355;
    local($lr1)=$_[2]/2.355;
    local($lr)=$_[3]/2;
    local($lA)=(1-$lSR)*$lr0*$lr0/$lr1/$lr1;
# Encircled Energy Ratio= (SR ∫_0^r exp(-r^2/r_0^2) dr^2 + A ∫_0^r exp(-r^2/r_1^2) dr^2 )/∫exp(-r^2/r_0^2) dr^2 
#                       = SR (1-exp(-r^2/r_0^2))+ A*r_1^2/r_0^2*(1-exp(-r^2/r_1^2))
    return($lSR*(1-exp(-$lr*$lr/$lr0/$lr0)) + $lA * $lr1*$lr1/$lr0/$lr0*(1-exp(-$lr*$lr/$lr1/$lr1)));
}