#!/usr/local/bin/perl -w

# flux scale in mJy
%f0
  = (
     4000 => 4.6,
     4500 => 4.3,
     5000 => 3.9,
     5500 => 3.5,
     6000 => 3.2,
     6500 => 3.0,
     8000 => 2.5
    );

# sky brightness in uJy/arcsec^2
%f0_sky
  = (
     4000 =>  6.7,
     4500 =>  4.3,
     5000 =>  5.6,
     5500 =>  7.3,
     6000 => 10.6,
     6500 => 15.6,
     8000 => 20.1
    );

%atmospheric_transmittance
  = (
     4000 => 0.66,
     4500 => 0.79,
     5000 => 0.83,
     5500 => 0.85,
     6000 => 0.86,
     6500 => 0.90,
     8000 => 0.92
    );

%sky_brightness
  = (
     'Dark' , {
	       4000 => 22.6,
	       4500 => 22.4,
	       5000 => 22.0,
	       5500 => 21.5,
	       6000 => 20.7,
	       6500 => 19.8,
               8000 => 19.3
	      },
     'Grey' , {
	       4000 => 21.3,
	       4500 => 21.3,
	       5000 => 21.2,
	       5500 => 20.8,
	       6000 => 20.3,
	       6500 => 19.5,
               8000 => 19.2
	      },
     'Bright' , {
	       4000 => 19.0,
	       4500 => 18.8,
	       5000 => 18.6,
	       5500 => 18.5,
	       6000 => 18.5,
	       6500 => 18.7,
               8000 => 18.2 
	      }
    );

%read_noise
  = (
     'EEV' , {
              'XTRASLOW' => 2.7,
              'SLOW'     => 2.7,
              'NORMAL'   => 3.4,
              'FAST'     => 4.4
             }
    );

%ccd_efficiency = (
     'EEV' ,  {
                4000 => 0.732,
                4500 => 0.83,
                5000 => 0.917,
                5500 => 0.90,
                6000 => 0.88,
                6500 => 0.862,
                8000 => 0.54
               });

%dispersion
  = (
	       '1200B' => 0.80,
	       '1200V' => 0.80,
	       '1200R' => 0.80,
	       '1200I' => 0.80,
	       '600U'  => 1.56,
	       '600V'  => 1.56,
	       '600R'  => 1.56,
	       '300B'  => 3.05,
	       '316R'  => 2.90,
	       '270R'  => 3.39,
    );

%grating_efficiency
  = (
     '1200B' , {
		4000 => 0.580,
		4500 => 0.548,
		5000 => 0.518,
		5500 => 0.455,
		6000 => 0.398,
		6500 => 0.347,
                8000 => 0.20
	      },
     '1200V' , {
		4000 => 0.482,
		4500 => 0.574,
		5000 => 0.613,
		5500 => 0.622,
		6000 => 0.601,
		6500 => 0.559,
                8000 => 0.448
	      },
     '1200R' , {
		4000 => 0.100,
		4500 => 0.220,
		5000 => 0.274,
		5500 => 0.410,
		6000 => 0.494,
		6500 => 0.545,
                8000 => 0.58
	      },
     '1200I' , {
		4000 => 0.000,
		4500 => 0.000,
		5000 => 0.119,
		5500 => 0.144,
		6000 => 0.191,
		6500 => 0.216,
                8000 => 0.36
	      },
     '600U' , {
		4000 => 0.497,
		4500 => 0.418,
		5000 => 0.362,
		5500 => 0.314,
		6000 => 0.256,
		6500 => 0.000,
                8000 => 0.00
	      },
     '600V' , {
		4000 => 0.613,
		4500 => 0.694,
		5000 => 0.738,
		5500 => 0.726,
		6000 => 0.682,
		6500 => 0.624,
                8000 => 0.43
	      },
     '600R' , {
		4000 => 0.000,
		4500 => 0.000,
		5000 => 0.467,
		5500 => 0.560,
		6000 => 0.636,
		6500 => 0.641,
                8000 => 0.595
	      },
     '300B' , {
		4000 => 0.73,
		4500 => 0.70,
		5000 => 0.67,
		5500 => 0.59,
		6000 => 0.51,
		6500 => 0.45,
                8000 => 0.30
	      },
     '316R'  , {
                4000 => 0.35,
                4500 => 0.56,
                5000 => 0.70,
                5500 => 0.70,
                6000 => 0.73,
                6500 =>	0.68,
                8000 => 0.54 
              },
     '270R' , {
		4000 => 0.300,
		4500 => 0.426,
		5000 => 0.596,
		5500 => 0.670,
		6000 => 0.748,
		6500 => 0.610,
                8000 => 0.645
	      }

    );


# test output

# ($snr,$nobj,$nsky) = 
#   signal_to_noise('EEV',
# 		  'SLOW',
# 		  '300B', 
# 		  30 ,
# 		  'Dark',
# 		  16 ,
# 		  4500 ,
# 		  1);

# print "snr is $snr \n";
# print "sky is $nsky \n";
# print "nobj is $nobj \n ";

#
# now do the actual work
#
# calculate the object and sky flux densities and the atmospheric
# transmittance

sub signal_to_noise {
  (
   $ccd,
   $speed,
   $grating,
   $zd,
   $moon,
   $mag,
   $wavel,
   $exptime,
  ) = @_;
  
  my ($SNR,$nobj,$nsky,$n_dn);
  
  $lines = "$grating";                          # get lines/mm from grating
  $lines =~ s/[A-Z]//;                          # (used for vignetting)

  $fnu  = $f0{$wavel} * 10.**(-0.4*($mag-15.));
  
  $fsky = $f0_sky{$wavel}
  * 10.**(-0.4*($sky_brightness{$moon}{$wavel}
		- $sky_brightness{Dark}{$wavel}));

  $disp=$dispersion{$grating};

  # get area of seeing disk in slit
  $area = 0.95;
  # convert ZD to airmass, converting to radians first
  $secz = 1./cos($zd*0.01745329);
  
  # get efficiency: spectrograph * grating * CCD
  
  $tseff=0.15;

  $seff = $grating_efficiency{$grating}{$wavel}
  * $ccd_efficiency{$ccd}{$wavel};
  
  $eff = $tseff * $seff;

  # now calculate the number of photons
  $nobj = 1.5e4 * 1.00 * $fnu * $area / $wavel
    * $atmospheric_transmittance{$wavel}**$secz * $eff;

  $nsky = 1.5e4 * 1.00 * $fsky*1.e-3 * 28.3
    * $eff / $wavel ;
  
  $nobj = $nobj*$disp*$exptime;
  $nsky = $nsky*$disp*$exptime;
  #
  #now get the SNR
  $n_dn = $read_noise{$ccd}{$speed} * sqrt(4);
  $SNR = $nobj / sqrt( ($nobj + $nsky) + $n_dn**2 );

  return ($SNR,$nobj,$nsky,$n_dn,$eff);
}