/*
 * File: sn_d_t.c
 *
 * S/N (Direct Imaging) Table
 * - -  -               -
 *
 * Author: J.Brewer 1998
 * email: jbrewer@eso.org
 *
 * This program generates an HTML page with a table containing S/N
 * values for stars with a range of magnitude and a range of
 * exposure times.
 * 
 * Version Information:
 * 1.0 (April 1998): First release.
 *
 */

#include "lasilla.h"    /* La Silla Extinction coefficients */
#include "inst.h"	/* Instrumental parameters (defined by Makefile) */

#include <stdio.h>
#include <stdlib.h>	/* Needed for atof() function  */
#include <string.h>	/* Needed for strcmp() function */
#include <math.h>	/* Needed for pow() function */

#define NARGS 9		/* Number of program arguments */
#define NXT 50		/* Number of exposure times  */
#define PI 3.142	/* Self Explanatory... */


main(int argc, char *argv[])
{

  /* Arguments */
  float see, msky, l_mag, h_mag, airm, h_xt;
  int pnt_src;

  /* Indices and work values */
  float counts, extinc;
  float disk, skyflux, mag, magC, i, xt, objsignal, skysignal, noise;

  /* Were there enough arguments? */
  if (argc != NARGS)
    {
      char usage[]="filter seeing msky l_mag h_mag pnt_src airm h_xt";
      fprintf(stderr, "usage: %s %s\n", argv[0], usage);
      exit(-1);
    }

  /* Process command line arguments */
  *argv[1]= toupper(*argv[1]);	/* Filter, -> uppercase 	*/
  see	= atof(argv[2]);	/* Seeing, in arcseconds	*/
  msky	= atof(argv[3]);	/* Sky magnitude, in mag/sq"	*/
  l_mag	= atof(argv[4]);	/* Low magnitude for table	*/
  h_mag	= atof(argv[5]);	/* High magnitude for table	*/
  pnt_src	= atoi(argv[6]);	/* Pnt src (1) or ext obj (0)	*/
  airm    = atof(argv[7]);	/* Airmass (=1 at zenith)       */
  h_xt	= atof(argv[8]);	/* Max exposure time		*/

  /* Use filter to determine count rate and extinction */
  if	(!strcmp(argv[1], "U")){
    counts = U_CNT;
    extinc = U_EXT;
  }else if (!strcmp(argv[1], "B")){
    counts = B_CNT;
    extinc = B_EXT;
  }else if (!strcmp(argv[1], "V")){
    counts = V_CNT;
    extinc = V_EXT;
  }else if (!strcmp(argv[1], "R")){
    counts = R_CNT;
    extinc = R_EXT;
  }else if (!strcmp(argv[1], "I")){
    counts = I_CNT;
    extinc = I_EXT;
  }else{
    fprintf(stderr, "%s -- unknown filter\n", argv[1]);
    exit(1);
  }

  if(pnt_src){
    /* Calculate seeing disk in pixels */
    disk = PI*((see/SCALE)*(see/SCALE));
    /* Calculate count rate from sky, per pixel, per second */
    skyflux = disk*SCALE*SCALE*counts*(pow(10,(0.4*(15-msky))));
  }else{
    /* Calculate seeing disk in pixels (`seeing'=1 arcsec**2) */
    disk = (1/SCALE)*(1/SCALE);
    /* Calculate count rate from sky, per arcsec**2, per second */
    skyflux = counts*(pow(10,(0.4*(15-msky))));
  }

  /* Here's the HTML output... */
  /* Table with input parameters */
  printf("<!--#include file=\"T_head.html\"--> \n\n"); /* SSI directive */
  printf("<P> Object type: "); 
  if(pnt_src)
    printf("<B>point source.</B></P>\n\n");
  else
    printf("<B>extended object.</B></P>\n\n");
  printf("<TABLE BORDER=2>\n");
  printf("<TR><TH ALIGN=LEFT> Instrument    <TD> %s    \n", INST);
  printf("<TR><TH ALIGN=LEFT> CCD           <TD> %s    \n", CCD);
  printf("<TR><TH ALIGN=LEFT> RON (e-)      <TD> %4.2f \n", RON);
  printf("<TR><TH ALIGN=LEFT> ''/pixel      <TD> %4.2f \n", SCALE);
  printf("<TR><TH ALIGN=LEFT> Band          <TD> %s \n", argv[1]);
  printf("<TR><TH ALIGN=LEFT> e-/s for m=15 in %s filter \
<TD> %5.0f \n", argv[1], counts);
  printf("<TR><TH ALIGN=LEFT> Extinction in %s band \
<TD> %4.3f \n", argv[1], extinc);
  printf("<TR><TH ALIGN=LEFT> Seeing        <TD> %3.1f \n", see);
  printf("<TR><TH ALIGN=LEFT> Sky Magnitude <TD> %4.1f \n", msky);
  printf("<TR><TH ALIGN=LEFT> Min Magnitude <TD> %4.1f \n", l_mag);
  printf("<TR><TH ALIGN=LEFT> Max Magnitude <TD> %4.1f \n", h_mag);
  printf("<TR><TH ALIGN=LEFT> Airmass       <TD> %4.1f \n", airm);
  printf("<TR><TH ALIGN=LEFT> Max Exposure  <TD> %4.1f \n", h_xt);
  printf("</TABLE>\n\n");
  printf("<BR><BR>\n\n");
  /* Write out column headers for the main table */
  printf("<TABLE BORDER=1>\n");
  printf("<TH>        ");
  for(mag=l_mag; mag<=h_mag; printf("<TH>%-10.1f", mag++))
    ;
  /* Loop through exposures (new exposure for each row) */
  xt = 0;
  for(i=0; xt<=h_xt; ++i)
    {
      /* Start a new row, and print out xt header */
      printf("\n<TR>");
      xt = pow(2.,(i/2.));
      if(RON > (sqrt(skyflux*xt)))
	printf("\n<TH>*%-8.1f", xt);	/* RON dominated */
      else
	printf("\n<TH>%-8.1f", xt);	/* Sky dominated */
      /* Calculate S/N along the row */
      for(mag=l_mag; mag<=h_mag; ++mag)
	{
	  /* Correct for atmosperic extinction */
	  magC = mag + (airm-1)*extinc;

	  /* Calculate count rate from star, per second */
	  objsignal = counts*(pow(10,(0.4*(15-magC))))*xt;
	  skysignal = skyflux*xt;
	  noise     = sqrt(objsignal+skysignal+disk*(RON*RON));
	  printf("<TD>%-10.1f", (objsignal/noise));
	}
    }
  printf("\n</TABLE>\n\n");
  printf("<!--#include file=\"T_tail.html\"--> \n"); /* SSI directive */

  exit(0);
}