/* * File: sn_d_g.c * * S/N (Direct Imaging) Graph * - - - - * * Author: J.Brewer 1998 * email: jbrewer@eso.org * * This program will use the sm library to a plot of S/N against * exposure time. Note that both a postscript and a GIF file will * be produced. * * Version Information: * 1.0 (April 1998): First release. * */ #include "lasilla.h" /* La Silla Extinction coefficients */ #include "inst.h" /* Dutch instrumental parameters */ #include #include /* Needed for atof() function */ #include /* Needed for strcmp() function */ #include /* Needed for pow() function */ #define NARGS 9 /* Number of program argments */ #define STEPS 500 /* Number of points to calculate on graph */ #define XZERO 17000. /* Right of label cluster */ #define YZERO 8000. /* Bottom of label cluster */ #define PI 3.142 /* Self Explanatory... */ static float xt[STEPS], StoN[STEPS]; /* Used by subroutines */ static float counts, extinc; /* Used by subroutines */ static int pnt_src; /* Used by subroutines */ main(int argc, char *argv[]) { /* Functions Prototypes */ void do_plot(int gif, char *labels[]); /* Arguments */ float see, msky, mobj, airm, l_xt, h_xt; /* Indices and work values */ int i; float disk, skyflux, objflux, objsignal, skysignal, noise; /* Were there enough arguments? */ if (argc != NARGS) { char usage[]="filter seeing msky mobj pnt_src airm l_xt h_xt"; fprintf(stderr, "usage: %s %s\n", argv[0], usage); exit(-1); } /* Process the other 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" */ mobj = atof(argv[4]); /* Magnitude */ pnt_src = atof(argv[5]); /* Pnt src (1) or ext obj (0) */ airm = atoi(argv[6]); /* Airmass (=1 at zenith) */ l_xt = atof(argv[7]); /* Low xt for calculation */ h_xt = atof(argv[8]); /* High xt for calculation */ /* 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); } /* Correct mobj for atmospheric extinction */ mobj += (airm-1)*extinc; /* Calculate count rate from object, per second */ objflux = counts*(pow(10,(0.4*(15-mobj)))); if (pnt_src) { /* Calculate seeing disk in pixels */ disk = PI*((see/SCALE)*(see/SCALE)); /* Calculate count rate from sky in seeing disk, 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)))); } /* Calculate arrays holding xt and S/N values */ for(i=0; i < STEPS; i++) { xt[i] = l_xt + i*((h_xt-l_xt)/(STEPS-1)); objsignal = objflux*xt[i]; skysignal = skyflux*xt[i]; noise = sqrt(objsignal+skysignal+disk*(RON*RON)); StoN[i] = (log10(objsignal/noise)); } /* Use `do_plot' subroutine to make graphs */ do_plot(1, &argv[0]); /* Make the GIF file */ do_plot(0, &argv[0]); /* Make the PS file */ exit(0); } void do_plot(int gif, char *label[]) { char tmp[20]; float glimit(int horl, int count, float array[]); sm_defvar("TeX_strings","1"); if(gif) sm_device("GIF TEMP.gif"); else sm_device("postencap TEMP.ps"); sm_graphics(); sm_ctype("black"); sm_location(5000, 30000, 5000, 30000); sm_expand(1.2); sm_limits(glimit(0,STEPS,xt), glimit(1,STEPS,xt), glimit(0,STEPS,StoN), glimit(1,STEPS,StoN)); sm_ticksize(0., 0., -1., 10.); sm_box(1,2,0,0); sm_expand(1.8); sm_xlabel("{\\it Exposure Time (s)}"); sm_ylabel("{\\it S/N Ratio}"); sm_ctype("red"); sm_conn(xt, StoN, STEPS); /* Labels... */ sm_ctype("blue"); /* General labels... */ sm_expand(1.6); sm_grelocate(6000., 28000.); sm_putlabel(6, INST); sm_expand(0.9); /* User specified parameters... */ sm_grelocate(XZERO, YZERO+5000.); sm_putlabel(4,"Filter = "); sm_grelocate(XZERO, YZERO+5000.); sm_putlabel(6,label[1]); sprintf(tmp, "%-6.0f", counts); sm_grelocate(XZERO, YZERO+4000.); sm_putlabel(4,"e/s @ m=15 = "); sm_grelocate(XZERO, YZERO+4000.); sm_putlabel(6,tmp); sprintf(tmp, "%-6.3f", extinc); sm_grelocate(XZERO, YZERO+3000.); sm_putlabel(4,"Ext (mag/airm) = "); sm_grelocate(XZERO, YZERO+3000.); sm_putlabel(6,tmp); sm_grelocate(XZERO, YZERO+2000.); sm_putlabel(4,"Seeing ('') = "); sm_grelocate(XZERO, YZERO+2000.); sm_putlabel(6,label[2]); sm_grelocate(XZERO, YZERO+1000.); sm_putlabel(4,"Sky (Mag/arcsec^2) = "); sm_grelocate(XZERO, YZERO+1000.); sm_putlabel(6,label[3]); sm_grelocate(XZERO, YZERO+0000.); if(pnt_src) sm_putlabel(4,"Mag (pnt src) = "); else sm_putlabel(4,"Mag (arcsec^{-2}) = "); sm_grelocate(XZERO, YZERO+0000.); sm_putlabel(6,label[4]); sm_grelocate(XZERO, YZERO-1000.); sm_putlabel(4,"Airmass = "); sm_grelocate(XZERO, YZERO-1000.); sm_putlabel(6,label[6]); /* Telescope/Instrument parameters... */ sm_grelocate(XZERO+8500., YZERO+5000.); sm_putlabel(4,"CCD = "); sm_grelocate(XZERO+8500., YZERO+5000.); sm_putlabel(6,CCD); sprintf(tmp, "%-4.2f", RON); sm_grelocate(XZERO+8500., YZERO+4000.); sm_putlabel(4,"RON (e^-) = "); sm_grelocate(XZERO+8500., YZERO+4000.); sm_putlabel(6,tmp); sprintf(tmp, "%-4.2f", SCALE); sm_grelocate(XZERO+8500., YZERO+3000.); sm_putlabel(4,"Scale (''/pix) = "); sm_grelocate(XZERO+8500., YZERO+3000.); sm_putlabel(6,tmp); sprintf(tmp, "%-6d", FULLWELL); sm_grelocate(XZERO+8500., YZERO+2000.); sm_putlabel(4,"FullWell (ADU) = "); sm_grelocate(XZERO+8500., YZERO+2000.); sm_putlabel(6,tmp); sprintf(tmp, "%-6d", BIAS); sm_grelocate(XZERO+8500., YZERO+1000.); sm_putlabel(4,"BIAS (ADU) = "); sm_grelocate(XZERO+8500., YZERO+1000.); sm_putlabel(6,tmp); sprintf(tmp, "%-6.5f", CRRATE); sm_grelocate(XZERO+8500., YZERO+0000.); sm_putlabel(4,"CRs (/pix/hr) = "); sm_grelocate(XZERO+8500., YZERO+0000.); sm_putlabel(6,tmp); sm_grelocate(XZERO+8500., YZERO-1000.); sm_putlabel(4,"Params Updated: "); sm_grelocate(XZERO+8500., YZERO-1000.); sm_putlabel(6,UPDATE); sm_gflush(); /* Flush graphics */ sm_hardcopy(); /* Make a hardcopy */ } float glimit(int horl, int count, float array[]) { /* Return high (horl=1) or low (horl!=1) limit for a vector. Note that the limit is min/max value -/+ 5% of the range */ int i; float hi, lo; hi=lo=array[0]; for(i = 1; i< count; i++) { lo = (lo>array[i]) ? array[i]:lo; hi = (hi