#! /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)

#$BJ*M}Dj?t(B
$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
# $BBg5$$NF)2aN($,(B20%$B0J2<$NNN0h$O4QB,$G$-$J$$$H$9$k!#(B($thrputmisc[$i]=1e-50$B$K(B)$B!#(B
for( $i=0; $i<$#filter+1 ; $i++ ){
    $thrputmisc[$i]=1.0; # $B$^$:$O(Bthrputmisc$B$9$Y$F$,(B1$B$@$H$9$k!#(B
    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;}
}

#$BGX7JJ|<M$K$h$kC10L;~4VEv$j$K8!=P$5$l$k8w;R$N?t(B
for( $i=0; $i<$#filter+1 ; $i++ ){
#$B$^$:$O%U%i%C%/%9(B
    $f=$zmflux[$i]*(exp(log(10.0)*(-0.4*$bgmag[$i])));
#$B8!=P$5$l$k(BPhoton$B?t(B(/pix/sec)
    $ssky[$i]=$f*$diameter*$diameter/4*$PI*$dlambda[$i]*1e6*$thrput[$i]*$thrputinst[$i]/$H/$nu[$i]*$pixscale[$i]*$pixscale[$i];
}

#$BC10l@QJ,;~4V$O(B
for( $i=0; $i<$#filter+1 ; $i++ ){
    if($ssky[$i]==0) { printf("%d %f\n",$i,$ssky[$i]); }
#$B$^$:$O$I$N$/$i$$$G(Bwell$B$,$$$C$Q$$$K$J$k$+(B
#$B$b$7$b(Bsexpmin$B$@$H(Bwell$B$,0lGU$K$J$C$F$7$^$&$h$&$@$H!"$=$NGHD9$N4QB,$O$G$-$J$$$3$H$K$9$k(B($thrputmisc[$i]=1e-60$B$K(B)$B!#(B
    if($lambda[$i]<1e-6){
# $B2D;k$G$O(B5$BJ,(B
	$sexposure[$i]=300
    }
    elsif($lambda[$i]<2.5e-6){
# $B%-%s%;%-%,%$$G$O(B2$BJ,(B
	$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$B$,(B3600$BIC$r1[$($k$N$G$"$l$P6/@)E*$K(B3600$BIC$K$7$F$7$^$&!#(B
    if($sexposure[$i] > 3600) {
	$sexposure[$i]=3600;
    }

#$B$D$.$K!"$I$N$/$i$$$G(Bbackground limited$B$K$J$k$+(B
# 4 Readnoise = BGnoise $B$,J,4tE@$G!"$=$l$KE~C#$9$k@QJ,;~4V$r(Bminexposure$B$H$9$k(B

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

#minexposure$B$,(B20$BIC0J2<$G$"$l$P6/@)E*$K(B20$BIC$K$7$F$7$^$&!#(B
    if($minexposure[$i] < 20) {
	$minexposure[$i]=20;
    }
#    print("$minexposure[$i] ");

##minexposure $B$H(B sexposure $B$rHf3S$7$F!">.$5$$J}$r:NMQ$9$k!#(B
    if( $minexposure[$i] < $sexposure[$i] ){
	$sexposure[$i]=$minexposure[$i];
    }

#$BAm@QJ,;~4V$G(Bwell$B$,$"$U$l$J$1$l$P!D(B
    if($sexposure[$i] >= $exposure) {
	$sexposure[$i]=$exposure;
	$n[$i]=1;
    }
#$BAm@QJ,;~4V$G(Bwell$B$,$"$U$l$k$J$i!D(B
    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 $B$NJ,$rLa$9!#(B
    $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)$B$K$h$C$F!"$I$NDxEY$N8w$,$b$H$a$h$&$H$7$F$$$k%"%Q!<%A%c!<$KF~$C$F$/$k$N$+$r7W;;!#(B
# ($B40A4$J(Bdiffraction limit)=(diffraction core)+(seeing halo)
# $B"i(Bexp(-r^2/r_0^2) dr^2= SR*$B"i(Bexp(-r^2/r_0^2) dr^2 + A*$B"i(Bexp(-r^2/r_1^2)dr^2
#   $B&H(B_diffrac=2.3548 * r_0
#   $B&H(B_seeing=2.3548 * r_1
#   A = $B%7!<%$%s%0%O%m$N%T!<%/(B
# $B&H$H(BSR$B$,J,$+$C$F$$$F!"(BA$B$r5a$a!"$5$i$K$"$k%"%Q!<%A%c!<Fb$KF~$C$F$$$k%(%M%k%.!<$N3d9g$r;;=P!#(B
# $B$^$:$O(BA$B$r5a$a$k%5%V%k!<%A%s(B
sub halo_peak {
# $_[0]=strehl ratio
# $_[1]=r_0
# $_[2]=r_1
    local($lSR)=$_[0];
    local($lr0)=$_[1];
    local($lr1)=$_[2];
# $B"i(Bexp(-r^2/r_0^2) dr^2= 2$B&P"i(Bexp(-r^2/r_0^2) r dr=$B&P(B * r_0^2 $B$h$j(B
# A=(1-s)(r_0/r_1)^2
    return((1-$lSR)*$lr0*$lr0/$lr1/$lr1);
}

# $B<!$K(BSR,A,$B&H(B_0, $B&H(B_1$B$+$i!"&H$ND>7B$N(Bencircled energy ratio$B$r$b$H$a$k(B
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 $B"i(B_0^r exp(-r^2/r_0^2) dr^2 + A $B"i(B_0^r exp(-r^2/r_1^2) dr^2 )/$B"i(Bexp(-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)));
}

#$B>eFs$D$r$^$H$a$?$b$N(B
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 $B"i(B_0^r exp(-r^2/r_0^2) dr^2 + A $B"i(B_0^r exp(-r^2/r_1^2) dr^2 )/$B"i(Bexp(-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)));
}
