/* Simple program to compute the Signal to Noise ratio of a telescope + CCD setup

 

Author       :      Juan Lacruz, May 2003

References :      Detecting and measuring faint point sources with a CCD, Herbert Raab 2002.

How to calculate the share for the central pixel.htm

*/

 

#include <stdio.h>

#include <math.h>

 

main()

{

 

       float Pi = 3.1415926535;

       float MirrorDiameter;

       float FocalLenght;

       float FocalRatio;

       float Obstruction;

       float LightCollectingArea;

 

       float PixelSize;

       float PixelScale;

       float DarkCurrent;

       float ReadOutNoise;

       float QuantumEfficiency;

      

       float IntegrationTime;

       float FWHM;

       float Sigma;

 

       float x;

       float erfx;

 

       float ObjectMagnitude;

       float ObjectFlux;

       float ObjectSignal;

       float ShareForCentralPixel;

       float ObjectSignalInCentralPixel;

       float ObjectNoiseInCentralPixel;

       float SkyBackground;

       float SkyBackgroundFlux;

       float SkyBackgroundSignal;

       float SkyBackgroundNoise;

       float DarkSignal;

       float DarkNoise;

       float NoiseInCentralPixel;

       float PeakSNR;

      

 

       printf("\n    Signal to Noise Ratio calculation \n");

       printf("\n    Telescope + CCD camera system \n\n\n\n");

 

 

       printf("Mirror diameter (m): ");

       scanf("%f",&MirrorDiameter);

 

       printf("Focal lenght (m): ");

       scanf("%f",&FocalLenght);

 

       printf("Obstruction (m): ");

       scanf("%f",&Obstruction);

 

       printf("Pixel size (um): ");

       scanf("%f",&PixelSize);

 

       printf("Dark current (e-/sec/pixel): ");

       scanf("%f",&DarkCurrent);

 

       printf("Readout noise (e-/pixel): ");

       scanf("%f",&ReadOutNoise);

 

       printf("Quantum efficiency (%): ");

       scanf("%f",&QuantumEfficiency);

 

       printf("Integration time (s): ");

       scanf("%f",&IntegrationTime);

 

       printf("FWHM (arc sec): ");

       scanf("%f",&FWHM);

 

       printf("Object magnitude (m): ");

       scanf("%f",&ObjectMagnitude);

 

       printf("Sky background (mag / squared arc sec): ");

       scanf("%f",&SkyBackground);

 

 

 

 

       FocalRatio = FocalLenght/MirrorDiameter;

 

       printf("Focal ratio %f \n", FocalRatio);

 

       LightCollectingArea = Pi * (pow(MirrorDiameter/2,2) - pow(Obstruction/2,2));

 

       printf("Light collecting area %f (m^2)\n", LightCollectingArea);

 

       PixelScale=asin(PixelSize/1000000/FocalLenght)*360/2/Pi*3600;

      

       printf("Pixel scale %f arc sec/pixel\n", PixelScale);

 

 

       ObjectFlux = 4E10 /pow(2.5,ObjectMagnitude);

       ObjectFlux = ObjectFlux * IntegrationTime * LightCollectingArea;

 

       printf("Object flux %f photons\n", ObjectFlux);

      

       ObjectSignal = ObjectFlux * QuantumEfficiency/100;

       printf("Object signal %f (e-)\n", ObjectSignal);

 

       /* Now calculate share for central pixel */

 

       Sigma = FWHM/2.355;

       x = PixelScale/2/Sigma/sqrt(2);

       erfx = 2/sqrt(Pi)*(x - pow(x,3)/3+pow(x,5)/10-pow(x,7)/42);

       ShareForCentralPixel = pow(erfx,2);

 

       printf("Share for central pixel %f\n",ShareForCentralPixel);

 

       ObjectSignalInCentralPixel = ObjectSignal * ShareForCentralPixel;

 

       printf("Object signal in central pixel %f (e-)\n", ObjectSignalInCentralPixel);

 

       ObjectNoiseInCentralPixel = sqrt(ObjectSignalInCentralPixel);

 

       printf("Object noise in central pixel %f (e-)\n", ObjectNoiseInCentralPixel);

      

 

       SkyBackgroundFlux = 4E10/pow(2.5, SkyBackground);

       SkyBackgroundFlux = SkyBackgroundFlux * IntegrationTime * LightCollectingArea * pow(PixelScale,2);

      

       printf("Sky background flux %f (photons/pixel)\n", SkyBackgroundFlux);

 

       SkyBackgroundSignal = SkyBackgroundFlux * QuantumEfficiency / 100;

 

       printf("Sky background signal %f (e-/pixel)\n", SkyBackgroundSignal);

 

       SkyBackgroundNoise = sqrt(SkyBackgroundSignal);

 

       printf("Sky background noise %f (e-/pixel)\n", SkyBackgroundNoise);

 

       DarkSignal = DarkCurrent * IntegrationTime;

      

       printf("Dark signal %f (e-/pixel)\n", DarkSignal);

 

       DarkNoise = sqrt(DarkSignal);

 

       printf("Dark noise %f (e-)\n", DarkNoise);

 

       NoiseInCentralPixel = pow(ObjectNoiseInCentralPixel,2) + pow(SkyBackgroundNoise,2);

       NoiseInCentralPixel = NoiseInCentralPixel + pow(DarkNoise,2) + pow(ReadOutNoise,2);

       NoiseInCentralPixel = sqrt(NoiseInCentralPixel);

 

       printf("Noise in central pixel %f (e-)\n", NoiseInCentralPixel);

 

       PeakSNR = ObjectSignalInCentralPixel/NoiseInCentralPixel;

 

       printf("Peak Signal to noise ratio %f  \n\n\n\n", PeakSNR);

 

}