C ~crb/ing/signal.f C C Chris Benn La Palma crb@ing.iac.es C C Predicts photon counts from target and from sky, and resulting C signal-to-noise, for spectroscopic and direct-imaging observations. C This version is for the instruments of the Isaac Newton Group, La Palma, C but only the DATA statements are site- or instrument-specific. C C If you are viewing this code with a web browser, you will probably C find chunks of concatenated lines near the end of the program. C This is triggered by html 'pre' tags embedded in write statements (used C when SIGNAL is run via the web interface), but if the code is saved C to a file, it should look normal. C C COMPILING AND LINKING C --------------------- C To compile and link under unix, use: C f77 signal.for -o signal.exe C (option -e needed if modifications extend any lines beyond 80 characters). C The program is self-contained, no libraries or other files needed (all C data are included in the DATA statements). No graphics. C If you have any problems running this code on your machine, please let me C know. C The program is written in Fortran 77, with the following widely-used C extensions (as flagged by f77 -ansi): IMPLICIT NONE, DO...ENDDO, C subroutine names > 6 characters or including '_', C non-standard 'name=' (why?), non-standard comment (separating data C statement and data?). Code checked with ftnchek (for date of last C check, see history section of comments). C C If you obtained the executable by ftp, you may need C chmod +x signal C to reset the protection on the file, before running it. C C C WHAT THE PROGRAM DOES C --------------------- C The program calculates the number of photons/A from the object C (SOBJ) and photons/A/arcsec**2 from the sky (SSKY) as = C 10**(-mag/2.5) where mag is OBJM or SKYMAG(IBAND) C OBJM is magnitude for a point object, C mag/arcsec**2 for extended object C *TIME integration time in sec C *PHOTONS(IBAND)*10000 conversion mag -> C photons/m**2/sec/A C *EXTIN(IBAND)**AIRMASS zenith transmission by atmosphere C (implemented from v11.0) C *TELAREA(ITEL(IINSTR)) telescope area (m**2) C *EFFMIRROR**IMIRROR(IINSTR) efficiency of main mirrors C *VIGN fraction of light passed by slit, C fibre or lesnlet (spectroscopic obs C of point sources only) C *EFFINSTR(IBAND,IINSTR) throughput of instrument (including C derotator e.g Nasmyth and corrector C e.g. PF, if any, and excluding other C components listed here) C *EFFGRAT(IBAND,IGRAT,IINSTR) efficiency of grating C *EFFDET(IBAND,IDET) q. efficiency of detector C *TRUETHRU(IBAND,IINSTR) ratio empirical/theoretical throughput, C now determined for most ING instruments C (see SPECTRUM March 1995), values C typically 0.7 - 1.1, other values C extrapolated or set to 1.0 (all were C 1.0 for pre-1995 versions of SIGNAL). C C If IPX=0 (point-source calculation), OBJM is the apparent mag, C if IPX=1 (surface-brightness calc.), OBJM is apparent mag/arcsec**2 C ITYPE = 1 for an imaging calculation, 2 for spectroscopic, C 3 for integral-field or fibre spectroscopy (in which case IPX is C forced to be 0). C C (1) IMAGING - POINT SOURCE (ITYPE=1, IPX=0) C Number of photons from object (SOBJ2) and sky/pixel (SSKY2) are C calculated as: C SOBJ2=SOBJ*BAND(IBAND) C SSKY2=SSKY*BAND(IBAND)*ARCPIX*ARCPIX C And S:N ratio for an aperture of diameter 2 FWHM is: C SNR=SOBJ2/SQRT(SOBJ2+SEE*(SSKY2+READ**2)) C where SEE is pi*(FWHM/ARCPIX)**2 pixels (number of pixels over C which integrating) C C (2) IMAGING - EXTENDED SOURCE (ITYPE=1, IPX=1) C Number of photons from object/pixel (SOBJ2) and sky/pixel (SSKY2) are C calculated as: C SOBJ2=SOBJ*BAND(IBAND)*ARCPIX*ARCPIX C SSKY2=SSKY*BAND(IBAND)*ARCPIX*ARCPIX C And S:N ratio per pixel is: C SNR=SOBJ2/SQRT(SOBJ2+SSKY2+READ**2) C i.e. SEE is assumed to be 1 pixel C C (3) SPECTROSCOPY - POINT SOURCE (ITYPE=2, IPX=1) C Number of photons per pixel-step in wavelength (SOBJ2) and C number of photons/pixel from the sky (SSKY2) are calculated as: C SOBJ2=SOBJ*DISPPIX C SSKY2=SSKY*SLIT*ARCPIX*DISPPIX C And S:N ratio per pixel-step in wavelength (i.e. per pixel of C extracted spectrum) is: C SNR=SOBJ2/SQRT(SOBJ2+SEE*(SSKY2+READ**2)) C where SEE is assumed to be 2*(FWHM/ARCPIX) C C (4) SPECTROSCOPY - EXTENDED SOURCE (ITYPE=2, IPX=2) C Number of photons/pixel from object from the object (SOBJ2) and C photons/pixel from the sky are calculated as: C SOBJ2=SOBJ*SLIT*ARCPIX*DISPPIX C SSKY2=SSKY*SLIT*ARCPIX*DISPPIX C And S:N ratio per pixel is: C SNR=SOBJ2/SQRT(SOBJ2+SSKY2+READ**2) C i.e. SEE is assumed to be 1 pixel C C (5) INTEGRAL-FIELD / FIBRE SPECTROSCOPY (ITYPE=3, IPX=1) C Number of photons per pixel-step in wavelength (SOBJ2) and C number of photons/pixel from the sky (SSKY2) are calculated as: C SOBJ2=SOBJ*DISPPIX C SSKY2=SSKY*SLIT*SLIT*3.14159/4.*DISPPIX C (i.e. phot/pixel-step, not per pixel) C And S:N ratio per pixel-step in wavelength (i.e. per pixel of C extracted spectrum, for one lenslet or fibre) is: C SNR=SOBJ2/SQRT(SOBJ2+SSKY2+SEE*READ**2) C where SEE is assumed to be 3 pixels (spread of light orthogonal C to dispersion direction) C C The program also calculates the scales at the detector in arcsec/pixel C (ARCPIX) and A/pixel (DISPPIX) as: C ARCPIX=ARCMM(IINSTR) *PIXEL(IDET) C DISPPIX=DISP(IGRAT,IINSTR)*PIXEL(IDET) C C Where no empirical data are available, SIGNAL uses interpolated, C extrapolated or (as a last resort) theoretical throughputs. C Empirical throughputs are typically 80 - 90% of theoretical. C NB some of the instruments and detectors documented by SIGNAL are not C currently offered by ING (e.g. WHT FOS, JKT RBS, IPCS). C C The program loops, at each loop either obtaining the parameters C interactively from the user, via the menu, or reading them C non-interactively from the data file SIGNAL.INPUT (one loop only C if this file has been generated from the web interface). C C LIMITATIONS (some of this is summarised in subroutine helpsub) C ----------- C The program does not take into account the following: C - Variation of atmospheric (non-grey) and dust (grey) extinction. C The default values are from La Palma technical note 31. C La Palma suffers from Saharan dust in the summer, extinction typically C 0.2 - 0.4 mag, up to 1 mag. The dust is grey but the extinction C varies from night to night (see the ING web page for C a nightly record from the Carlsberg Meridian Circle). C - Degradation of the main telescope mirrors between aluminisings. This may C be up to 10%, more in the UV. C - Vignetting and optical aberrations (broadening the image) at large C field radius (see the ING Observers' Guide for more information). C - Vignetting by the spectrograph slit, important for point sources if C seeing ~ slit width, as is usually the case, and wavelength dependent C if slit position angle is not = parallactic angle. C The program light_in_slit (see ING web page) can be used to calculate C the magnitude of the effect. C - Loss of light to colour, ND and polarisation filters in the light path. C But SIGNAL assumes that order-sorting filters HAVE been inserted where C necessary e.g. to exclude second-order light in ISIS red arm! C - Loss of grating efficiency at large angles, or when the grating is C reversed or overfilled C - Conversion between recorded counts from CCD and number of detected C photons (i.e. the effects of gain = e-/ADU, bias level, dark current, C non-linearity of response, cosmic rays etc.) C Quoted CCD readout noise and e-/ADU are for the commonest readout speed C - Brightening of the sky due to: C - twilight (if sun < 18 deg below horizon); C - moonlight (if moon not below horizon); C - solar cycle (values used are for solar minimum, sky brightens up to C 0.4 mag at solar maximum); C - zodiacal light (sky is 0.4 mag brighter on ecliptic than at poles); C - Milky Way; C - airmass (sky is 0.2 mag brighter at airmass 1.5); C - light pollution - the measured zenith values used probably include no C significant contibution from light pollution, but sky is C significantly brighter when very dusty (> 0.3 mag extinction) C - Signal-to-noise for spectroscopy of unresolved objects is given C per pixel step in wavelength, which may not be the same as the C wavelength resolution (typically 2 pixels). C Notes on ING-specific instruments, filters and detectors: C - ISIS polarisation optics have throughput 0.9 approx, ISIS dichroics C about 0.85 C - For FOS, U,B are for second order, V,R,I are for first order C - The BVRI filters are assumed to be Harris. Note that the throughput C of the KPNO B filter is about 25% greater than Harris C The U filter is the 50-mm type. The 125-mm WHT PF U filter has a C transmission about a factor of two lower. C C C HISTORY C ------- C 1990 Nov 26 Version 1.0 in VAX Fortran C includes non-interactive option 'X' C 1991 Jun 11 Version 2.0 C 1992 Feb 8 Corrected erroneous I, Z CCD sensitivities C 1992 Mar 12 Added detector QE for Tek CCD (version 3.0) C 1992 Mar 25 Updated GEC/EEV CCD responses with PRJ's latest data C (version 3.1) C 1992 Apr 3 Added LDSS spectroscopic and imaging C Generalised code to be instrument-independent (version 4.0) C 1993 Oct 14 Updated UES efficiencies using Max Pettini's data C 1993 Dec 9 Added INT FORD and GEC7 QE C 1994 Sep 6 Version 6.0: added WHT prime (instr 19), removed JKT QUBES C and grism options; added WHT WYFFOS and hitchhiker; C updated readout noises and default CCD; C corrected IPCS2 pixel size to 10 microns; C added derotator eff 0.75 for UES; C added R1200B grating (assume eff as R600B) C 1995 Jan 23 Added TRUETHRU empirical/predicted data (+ extrapolations) C WHT FOS throughput based on 1/1/92 measurement C Version 7.0 (-> UK/NL/ES site managers 6/3/95) C 1995 Mar 24 Removed incorrect warning that no throughput data for WYFFOS C Added correct DISP for WYFFOS C 1995 May 10 Removed tabs. Version 7.1 C 1995 May 16 Corrected WHT PF TRUETHRU(B) from 0.4 to 0.6. Version 7.2 C 1995 Nov 14 Revised (slightly upwards) aux-port TRUETHRU on basis of C throughput measurements of night 13/11/95. Version 7.3 C 1995 Nov 24 Added table showing status of empirical data C 1996 Aug 14 Broke up data statements to avoid more than 19 C continuation lines (Fortran 77 standard, adhered to by C some UNIX machines) C Version 8.0 C 1996 Sep 20 Removed carriagecontrol='list' which makes Linux choke C Added Loral QE C Version 9.0 (emailed to ING/RGO webmans + UK/NL/ES C site managers) C 1996 Oct 1 LORAL readnoise=8, version 9.1 C 1996 Oct 10 Added help comment about B(KP) vs B(Harris) C 1996 Nov 11 Updated sky-brightness values (version 9.2) C 1996 Dec 17 Added ISIS R grating efficiencies in Z (version 9.3) C 1997 Mar 3 Corrected WHT PF U effinstr to 1.0->0.6, truethru 0.7->1.1 C (version 9.4) C Mar 14 Added WYFFOS empirical throughput data from 11/96 comm run, C provisonal result that 0.5 mag poorer than ISIS (version 10) C Dec 2 Added EEV10 CCD C Removed dependence of ADU on INSTRUMENT C Fine-tuned sky-brightness values C 1998 Feb 13 Structural changes and enhancements not affecting C object and sky count rates for main instruments: C - Added calculations of S:N for spectrsocopy and for C extended images both imaging and spectroscopy C - Option to set effective bandwidth, useful if using C a filter with very different transmission from C default set (e.g. WHT PF U) or narrow band C - Added JHK bands, WHIRCAM InSb detector, WHIRCAM instrument C - Added STJ detector, human eye, scrapped IPCS and redundant C CCDs C - FWHM now refers to fwhm of star only, IPX takes over C function of flagging whether calculation is for star C or per arcsec**2 for extended object. C - Added calculation of detector sizes in both directions. C - Reorganised to make addition of new telescopes, C instruments etc. easier. Some variable names changed. C Extended arrays to allow 10 bands instead of 6 C and 50 instruments instead of 20. C - Made instrument-menu 2-column C - Help simplified, most additional info now in comments C in code C - Updated readout-noise and gain values, now reports both C - Corrected error in WHT FOS V throughput C - Added extinction-correction code, but set to 0.0 C - Made program check command issued by user C - Added option to change telescope area e.g. for JOSE, C Shack-H, which use sub-apertures C (version 10.2, checked ISIS + INT-imaging calcs give same C counts as for previous versions) C 1998 Feb 17 Changes which DO affect object and sky count rates: C - Implemented extinction correction C - Updated sky-brightness values (small changes) C - Updated effective filter bandwidths BAND(I) for Harris C filters and new U filters (still waiting for Harris I scan) C NB big change in U bandwidth, prob. offset by inclusion of C atmospheric ext. in U, so TRUETHRU not much different C for imaging. C - Changes to theoretical throughput values: C - UES throughput changed 0.36 to 0.45 (new derotator) C ISIS throughput changed 0.45 to 0.51 or 0.42 (two arms) C WYFFOS throughput changed 0.36 to 0.51 (error) C WHT PF B throughput changed 1.0 to 0.8 (PF corrector) C IDS500 throughput changed 0.41 to 0.44 (not same as 235) C - Multiplied INT PF U, B th. throughputs by 0.6, 0.8 because C of corrector. C - Updated table of empirical throughputs and adjusted C TRUETHRU values accordingly. C - Changed area over which sky flux counted from pi/4 FWHM**2 C to pi FWHM**2 for imaging calculations (Hans Deeg IAC C pointed out an inconsistency in previous versions, gave C correct counts from star and from sky, but an overoptimistic C estimate of resulting S:N for imaging, by integrating C sky over aperture of diameter FWHM, which doesn't include C all the flux) C (version 11.0) C 1998 Feb 20 - Changed to display band with wavelength when C a spectroscopic instr. (version 11.1) C 1998 Feb 21 - Changed to accept Oke AB mags too, altered some of C TRUETHRU to reflect difference between apparent mag C and Oke AB mag (version 11.2) C 1998 Apr 14 - Calculated object counts, for imaging of extended source, C corrected to be per pixel instead of per arcsec^2 C (error noted by Nigel Metcalfe) C - Two coding errors, accepted by SUN f77, rejected by C Richard McMahon's gnu g77 compiler: C (1) if(itype(iinstr.eq.1)) print *,... corrected to C if(itype(iinstr).eq.1) print *,... C (2) lines ending in comma, followed by lines starting C with comma, modified so only one comma C (even f77 -ansi doesn't catch this) C - IROT data statement accidentally commented out, C 90-deg rotation of EEV10a not being used C (error noted by Sally Oey) C - Minor changes to help + output format C (version 11.3) C 1998 Apr 20 Cleared up a few minor inconsistencies flagged by ftnchek C running on lpss20. Version 11.4 C 1998 Jun 3 Added CIRSI (INT or WHT). Version 11.5 C 1998 Jul 3 Added default detector numbers (IDEFDET) for each C instrument. Added INGRID. Version 11.6. C 1998 Jul 19 Tweaked sky-brightness values to match those given in C LP tech note 115 (Benn + Ellison). Version 11.7. C 1998 Aug 5 Adding sky-brightness calculator - estimates sky brightness C in each band with or without moonlight. Version 11.8. C Checked SIGNAL's output for ISIS and INT WFC agree previous C versions. Checked with ftnchek. C 1998 Sep 9 Adding lines at start of code to allow it to be used C in one of two modes after minor editing: C IP = 6 - usual interactive mode, writes to screen (unit = 6) C IP = 7 - reads parameters from signal.input, C writes to signal.output (unit=7), C which allows it be called from a web browser, C via the cgi (perl) script written by C Ashley James (summer student 1998) C Moved all parameter-setting code into two subroutines, C one for interactive, one for non-interactive. C Non interactive mode can also be triggered from menu, with C undocumented option X (untested). C Replaced 'print *' common to both modes by writes C (variable put free format in BUFFA, BUFFB) C Ruthlessly eliminated all data type logical C Version 12.0. C 1998 Sep 16 Minor changes to some format statements (error picked C up by Paul O'Brien's compiler). FTNCHEK OK. Version 12.1. C 1999 Mar 2 Updating data on Loral CCD so that it refers to the one in C INT FOS rather than to the defunct WHT ones. Version 12.2. C 1999 Mar 29 Correcting large-EEV read noise 11 -> 4 e-. Version 12.3. C 2000 Jun 21 Added some IDS gratings for WYFFOS. Version 12.4. C 2000 Sep 4 Removed print * formatting of output values from subroutine C free. C Changed K SKYMAG (i.e. sky + thermal) to 12.5 (for INGRID), C corrected INGRID gain to 5.6 e/ADU, QE changed from C estimate of 0.6 to .4, .55, .55, corrected ZJHK bandpasses C to be for INGRID. C Version 12.5 C 2000 Sep 13 New INGRID and UES throughput measurements, version 12.6. C 2001 Feb 12 New INGRID Ks skybr 11.8 mag/arcsec^2 from Johan Knapen's C measurements 4-12/00. C Rearranged instrument menu. C Added NAOMI/INGRID. C Added comment line to appear after o/p table for each instr. C Version 12.7 C 2001 Jul 19 New INGRID gain 5.3 (was 5.5), readnoise 32 (was 10), c scale 0.238 "/pixel, C J,H,K skybr values (fainter) and zeropoints (factor 1.6 C fainter, largely due to factor 1.4 increased throughput C for the new aspheric) C Version 12.8 C 2001 Oct 19 Commented K bgr high for NAOMI/INGRID C Version 12.9 C 2002 Feb 5 Added ITYPE(IINSTR) = 3 option for ifu/fibre spectroscopy, C so that now get sky counts and S:N for e.g. AF2/WYFFOS. C Added OASIS and MIT/LL CCD. C Added variable DEFSLIT - default slit width (or fibre C diam etc) for given instruments. C Corrected mismatches between instr names in .pl script C and html (including removal of multiple spaces, which C are not treated in a straightforward way). C Option to allow setting sky-br as D G or B. C Now outputs sky photons/A/sq.arcsec as well (McMahon's C suggestion 10/98). C Updated wyffos copy of data for ISIS gratings C (redward extrapolation had been done only for ISIS). C Version 13.0 C 2002 Feb 7 Added correction for vignetting by slit, fibre or C lenslet (VIGN). C Corrected parameter-order errors in sky-brightness C .pl script (noted by Ian Skillen 6/01). C Converted CBAND* and QBAND to 2-character, for future C implementation of additional filters. C Added variable SPREAD - number of pixels over which light C distributed orthognal to dispersion direction, for C fibre and lenslet spectrographs (2 for wyffos, 3 for C oasis). C Version 13.1 C 2002 Mar 2 Small changes to default dark, grey and bright sky C intensities. C Version 13.2 C 2003 Jun 4 Added NOT ALFOSC imaging + spectroscopy test C 2003 Jun 5 Changed number of detectors from 10 to 30, version 13.4 C Removed detector options: INT FOS Loral, old GEC C Added: Marconi 2148*4700, NOT CCD7 (Loral) C 2003 Jun 10 Adding NOT NOTCAM, MOSCA, version 13.5 C Removed superfluous print of readnoise to output text C Changed lengths and sizes of some arrays for C greater flexibility C Version 13.5 C 2003 Aug 25 Robert Greimel added mode ip=8, for use with his C new .php web interface. C Version 13.6 C 2013 Dec 18 Thanassis Katsiyannis. Adopted for ARISTARCHOS C Version AR1.0 C Old values of numerical parameters are often preserved as comment C lines (in the data block) C C C VARIABLES C --------- C ADU electrons per recorded count (DETECTOR) C AIRMASS airmass C ANS 1-char ans C ARCMM scale "/mm at detector for each INSTR C ARCPIX scale "/pixel at detector C AREA geometrical collecting area * mirror, instrument, grating C and detector efficiencies C BAND effective bandwidths of filters i.e. integral over T(l)dl C where T(l) is the transmission and l is wavelength = C area under filter transmission curve C BUFFA, BUFFB, BUFFC C 15-char buffers to which reals are written free format, C prior to output C CBAND, CINSTR, CGRATING, CDET - names of bands, C instr. combinations, gratings, detectors C CGRATING is (IGRAT,IINSTR) C CINSTRB buffer to hold instr menu before writing C DEFSLIT default slit width for each instrument C DISPMM dispersion A/mm for each grating (IGRAT,IINSTR) C DISPPIX dispersion A/pixel at detector C DX, DY sizes of detector in X and Y in " or A (whichever is C appropriate) C EXTIN extinction value in each band at 1 airmass C LOOP main loop counter C EFFDET detector efficiencies for each band (BAND,DETECTOR) C For GEC, sensitivities are for coated CCD C EFFGRAT grating efficiencies for each band (BAND,GRATING,INSTRUMENT) C EFFINSTR throughput of instrument for each band (BAND,INSTRUMENT) C EFFMIRROR reflective efficiency assumed for mirror C FLUXJY Jy flux (1E-26 W/Hz/m**2) for mag=0 in each band (Z is approx) C FWHM FWHM arcsec (seeing + object) C IAXIS normal direction on CCD of spectral dispersion (1 = x, 2 = y) C ('normal' is for TEK CCDs, EEV10a etc are 90-deg rotated) C IBAND, IINSTR, IGRAT, IDET - band, instrument, grating and detector C numbers C IDEFDET default detector number for each instrument C IDX,IDY axis labels: 0 = arcsec, 1 = A C IMIRROR number of mirrors between sky and instrument C IP 6 => interactive (writes to unit 6) C 7 => non-interactive reads from SIGNAL.DAT, writes to C SIGNAL.OUT (unit 7) C IPX calculation for stellar object (0) or per sq arcsec for C extended image (1) C IROT rotation of each detector: 0 means spectral direction C would be in X if mounted on ISIS, 1 means 90 deg rotated C ITEL telescope number for the given IINSTR (to get tel geom area) C 1 = JKT 2 = INT 3 = WHT 4 = EYE 5 = NOT C ITRANS determines order in which instrument menu written C ITYPE observation type for the given IINSTR C 1 = imaging, 2 = spectroscopic, 3 = ifu/fibre spectroscopic C MODE mode of use after parameters obtained: C 0 = S:N calc; 1 = sky br calc; 2 = help; 3 = skybr help; C 9 = end C OBJM magnitude of object C (per sq arcsec for spectroscopy with IPX=1) C NBAND, NINSTR, NGRATING(IINSTR), NDET - numbers of bands, C instruments, gratings and detectors C NPIXELS number of active pixels on CCD in direction (DETECTOR,AXIS) C PHOTONS photons/s/sq.cm/A for mag = 0 for each band C PIXEL pixel size (mm) for each detector C READ readout noise (electrons) C SEE 3.14159 * (FWHM/PIXSIZE)**2 i.e. number of pixels over which C 2D image is assumed to be spread (was 3.14159/4. * ... before C version 11.0), affects S:N ratio only C SKYMAG sky brightness in mag/arcsec**2 C SKYMAG0 set to SKYMAG at start, but SKYMAG0 cannot be changed C by the user (they are passed to subr SKYBRCALC) C SLIT, SLIT2 - slit-width (arcsec) C SNR signal-to-noise for imaging calculation C SOBJ number of photons/A from object C (per sq arcsec for spectroscopy with IPX=1) C SOBJ2 For spectroscopy, photons/pixel (sky or ext object with IPX=1) C or photons/pixel step in wavelength i.e. after integrating C in the spatial direction (IPX = 0) C SPREAD For fibre and lenslet spectrographs, the number of pixels C over which the spectrum is spread othrogonal to the C dispersion direction (typically 2 - 3). Note that this C assumes always same detector. C SSKY number of photons/arcsec**2/A from sky C SSKY2 SSKY per pixel for imaging or for normal spectroscopy C SSKY per pixel step in wavelength for IFU/fibre spec. C TELAREA geometric collecting area of telescope ITEL(IINSTR) C TIME integration time (sec) C TRUETHRU empirical/predicted throughputs (BAND,INSTR) C TRUE0 value of TRUETHRU passed to DISPLAY C UNIT(2) " or A (for writing out pixels scales in the two directions) C VERSION Version number and date C WAVE Wavelengths in microns of each band (used in v11.0 and after, C when calculating PHOTONS from FLUXJY C C Variables for passing to DISPLAY: C QINST,QBAND,QGRAT,QGRAT2,QDET,EXT, C GEOMAREA,THROUGHPUT,PIXSIZE,DETEFF, C PHOTONS2,READOUT,BWIDTH,SKY C C Variables used by sky-brightness calculator SKYBRCALC only: C SOL,ECLAT,GLAT,AIR,IM,ALPHA,R,ZM,AK C C ADDING NEW INSTRUMENTS, DETECTORS ETC. C -------------------------------------- C Data for new telescopes, instruments, gratings, detectors and bands C can be added by modifying the DATA statements (only). C In the table below are listed the principal data arrays, with C the array dimensions currently being: C B (bands) = 10 C G (gratings) = 10 C I (instruments) = 50 C D (detectors) = 30 C (e.g. 'G I' = one number required for each grating for each instrument). C C Variable Indices Meaning (see full list above for more detail) C C ADU D gain e-/count C ARCMM I "/mm C BAND B wavelength band C CBAND B band name C CDET D detector name C CGRATING G I grating name C CINSTR I instrument name C COMMENT I 1-line comment for instrument (printed by signal) C DEFSLIT I default slit width arcsec C DISPMM G I A/mm C EFFDET B D detector QE C EFFGRAT B G I grating throughput C EFFINSTR B I instrument throughput C EXTIN B extinction C FLUXJY B Jy for mag = 0 C IAXIS I 1 for dispersion along x, 2 for y C IDEFDET I default detector number C IMIRROR I number of telescope mirrors between instr and sky C IROT D (for ING instr only) C ITEL I telescope number C ITRANS I order in which instrument menu written C ITYPE I 1 = imager, 2 = spec, 3 = fibre or IFU spec C NGRATING I number of gratings C NPIXELS D no. pixels in x and y on detector C PIXEL D pixel size on detector C READ D read noise (e-) C SKYMAG B sky brightness C SPREAD I number of pixels over which spectrum spread perp C to wavelength, for fibre spec C TRUETHRU G I fudge factor to make prediction agree with C observed zeropoint C WAVE B wavelength in microns C C NBANDS, NGRATINGS, NINSTR and NDET may need changing. C If the detector arrays are changed, IDEFDET will need changing. C To add data involving a new telescope, may need to modify some or all C of IMIRROR, ITEL, TELAREA. C ITRANS will probably need reordering. C implicit none integer i,j,k,ninstr,ndet,iband,iinstr,igrat,ip,iin,iout,mode integer itrans(10),itel(10),itype(10),imirror(10),ngrating(10) integer iaxis(10),npixels(10,2),idefdet(10) integer idet,nband,loop,ipx,idx,idy real defslit(10),spread(10),vign real effinstr(10,10),effgrat(10,20,10),effdet(10,10) real truethru(10,10),dispmm(20,10),arcmm(10),telarea(10) real objm,time,snr,slit,arcpix,disppix,effmirror,area,sobj,ssky real sky,sobj2,ssky2,see,slit2,true0,fwhm,geomarea,ext,airmass real throughput,pixsize,deteff,photons2,bwidth,readout,dx,dy real pixel(10),read(10),adu(10),irot(10) real band(10),fluxjy(10),photons(10),skymag(10),skymag0(10), 1 extin(10),wave(10) character*1 ans character*2 cband(10),cbandlower(10),qband character*1 unit(2) character*6 qdet character*10 qgrat,qgrat2 character*15 buffa,buffb,buffc character*20 cgrating(20,10) character*25 cinstr(10),qinst character*28 version character*40 cdet(10) character*80 comment(50) C C Variables used by sky-brightness calaculator (only) C integer im real sol,eclat,glat,air,alpha,r,zm,ak data unit/'"','A'/ C C Version number and date C data version/'Version AR1.0 18 Dec 2013'/ C C Central wavelengths of bands used (from output of ~crb/ing/filter.f). C Typical wavelength ranges are given below. Actual values of effective C bandwidth (integral over T(l)dl, where T is transmission, l is wavelength) C are given in array BAND C BAND values are for: C U 3600 A (3300 - 3800 A) 50-mm square and 139-mm liq-CuSO4 U C B 4300 A (4000 - 5000 A) Harris C V 5500 A (5000 - 6000 A) Harris C R 6500 A (6000 - 7200 A) Harris C I 8200 A (7200 - 9200 A) Harris C Z 9500 A (8500 - 9700 A) RGO glass C J 12500 A (1.10 - 1.40 mic) WHIRCAM C H 16700 A (1.50 - 1.82 mic) WHIRCAM C K 21600 A (2.00 - 2.30 mic) WHIRCAM Kshort C data wave/3700.,4300.,5500.,6500.,8200.,0.,0.,0.,0.,0./ C data wave/3700.,4300.,5300.,5600.,7800.,0.,0.,0.,0.,0./ C data wave/3600.,4300.,5500.,6500.,8200.,9500., C 1 12500.,16700.,21600.,0.0/ C C Effective bandwidths UBVRI taken from INT prime-focus guide for C versions up to 10.2, and by integrating over UBVR filter curves with C ~crb/ing/filter.f for version 11.0 onward C NB still waiting for reliable scans of the I and Z filters C data band/750.,850.,950.,1500.,1650.,0.,0.,0.,0.,0./ C data band/1000.,2000.,2200.,3500.,4000.,0.,0.,0.,0.,0./! 1st estim. ARIST C 1 / 500., 720., 860.,1330.,1400., 700.,1600.,2900. ,3200. ,0./ ! >= 12.4 C C Mag -> Jy conversions from Bessell (1979 PASP 91 589) C and Bessell & Brett (1988 PASP 100 1134) C (prior to v11.0, were from Johnson H.L., 1966, ARAA, 4, 193) C data fluxjy C 1 /1880.,4440.,3810.,3010.,2430.,2200.,1570.,1020. , 640. ,0./ ! < 10.2 1 /1810.,4260.,3640.,3080.,2550.,2200.,1570.,1020. , 640. ,0./ ! >=11.0 C C C Photon rate for mag=0 (not used from v11.0 onward, since it can be C calculated from FLUXJY) taken from INT prime-focus manual C data photons C 1 / 783.,1463.,1050., 873., 376., 329., 194., 94.4, 44.1,0./ C New values are similar +-10% except R now 701, I now 449) C C Sky-brightness values are median dark-of-moon ING 1994-6 (solar C minimum) from Sara Ellison's summer 1996 work (for UBVRI; C Z is an interpolation between these and the known JHK values), C as reported in LP technical note 115. C K is for INGRID (Thomas Augusteijn email 7/8/00, previous values C were for WHIRCAM) C Sky brightness is similar (+- 0.2 mag) at other dark sites. C data skymag/21.7,22.3,21.3,20.6,18.9,0.0,0.0,0.0,0.0,0.0/ C 22.6,23.3,23.6,23.6,23.1,0.0,0.0,0.0,0.0,0.0/ C C Zenith extinction for La Palma, taken from the ING Observers' Guide C = constant (wavelength-dependent) gas + variable (grey) dust (= 0.05) C Used for version 11.0 onwards. C The values for other sites are similar. For CTIO, Frogel 1998, C PASP 110 200 gives .11, .06 and .09 in J, H and K. C data extin/0.62,0.25,0.12,0.07,0.04,0.0,0.0,0.0,0.0,0.0/ C data extin/0.62,0.17,0.09,0.04,0.01,0.0,0.0,0.0,0.0,0.0/ C C############################################################################## C Site/tel/instr-specific data start here and end at next line of hashes C C The TRUETHRU(BAND,INSTR) values are derived from the instrumental C zeropoints (apparent mag for 1 detected photon/sec for imaging C or Oke AB mag for 1 photon/sec/A for spectroscopy) given below. C Where two sets of values given, newest is used. C These measurements are likely to be C affected at the 10% level by changes in the primary-mirror C reflectivity between aluminisings. C Unless otherwise stated, the imaging data are through the 50-mm C CuSO4/glass U filters and the BVRI Harris set. C C Z E R O P O I N T S C Instrument U B V R I Date Source of information C ---------- ---- ---- ---- ---- ---- ---- --------------------- C WHT PF imaging 24.1 26.2 26.4 26.5 26.0 9409 Carter TEK C WHT aux-port 26.2 26.1 26.3 25.4 9511 Benn TEK C WHT TAURUS imag. 24.5 25.1 25.5 24.9 9203 GEMINI f/2, GEC C INT PF imaging 23.8 25.5 25.5 25.5 25.0 97? WEB Loral no. 4 C JKT imaging 21.0 23.1 23.3 23.5 23.0 9407 Peletier C (checked 21.1 23.1 23.2 23.4 23.0 9707 Benn TEK4) C NOT ALFOSC 23.6 25.7 25.7 25.5 24.9 0306 NOT page Loral C NOT MOSCA 24.6 26.2 26.0 25.7 25.2 0306 Telting Loral mosaic C NOTCAM HR + WF 0306 Telting Rockwell C J 23.8, H 23.8, K 23.3 C C WHT ISIS 18.1 18.2 17.6 18.5 17.4 9708 Smartt B/158B/LOR UBV C R/158R/TEK RI C 18.0 17.2 0101 Garcia B/158B/EEV12 C 18.0~18.0 R/158R/TEK4 C (no dichroic) C WHT FOS 17.0 18.1 17.2 9203 GEMINI C WHT LDSS spec 15.9 18.4 18.6 18.2 9404 Carter U (high), C B, V (blue), C R (red), TEK C WHT UES 16.3 17.0 17.0 16.7 15.5 9311 Pettini old derot, TEK C 16.7 16.9 16.9 16.0 0009 Telting new derot, SIT1 C 17.0 17.4 17.3 17.0 0101 Telting new derot, 2EEV C WHT WYFFOS refl. 16.5 16.5 16.5 9610? Bridges lower-limits C on throughput with C various gratings. C H2400B is factor 2 poorer. C 16.5 16.8 16.6 0107 Corradi small-fibres C R600B C WHT Hitchhiker 24.5 25.4 25.5 24.7 9206 GEMINI GEC C WHT WHIRCAM 9506 WHIRCAM users' guide C WHT CIRSI 9806 Richard McMahon C INT CIRSI 9806 Richard McMahon C WHT INGRID 0107 Almudena Zurita C 25.20 J, 25.28 H, 24.54 Ks C assuming 5.3 e-/count C INT IDS 235 13.8 16.3 16.4 16.2 15.1 97? WEB post-96-cleaning C R150V, TEK for BVRI C 9411 Harlaftis R150V TEK for U C 13.5 16.4 16.5 16.4 15.2 9810 Garcia EEV10 C 14.5 16.4 16.4 16.3 14.8 0010 Garcia TEK5 C INT FOS 13.7 15.7 15.7 17.0 15.9 9402 Harlaftis GEC C NOT ALFOSC 0306 Telting C C No empirical data for TAURUS FP, LDSS imaging, IDS500, C RBS or QUBES C C To obtain apparent mag, add the following numbers to an Oke AB mag: C C U B V R I C -0.8 0.2 0.0 -0.2 -0.4 C C The values of TRUETHRU below have been derived by judicious smoothing, C interpolation and extrapolation of the data above. Changes in the C instrument since the measurements were made (e.g. the new derotator C for UES have been taken into account. C In 3 cases, the derived TRUETHRU were significantly LARGER than 1.0. C ISIS R 1.7, ISIS I 1.4 are at risk of second-order contamination, C and have somewhat arbitrarily been set to 1.0. C The INGRID JHK values have been left as measured. C data ((truethru(i,j),i=1,10),j=1,10)/ C U B V R I Z J H K ? C 1 0.3,0.7,0.7,0.7,0.7,1.0,1.0,1.0,1.0,1.0, ! 10 IDS 235 1 0.7,0.7,1.0,1.0,1.0,1.1,1.0,1.0,1.0,1.0, ! 13 JKT imaging 1 0.0,1.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0, ! 22 Human eye 1 80*0.0/ data comment/24*' ', 1 ' (NB assumes read noise due to one IR-chip readout only)', ! 25 1 25 *' '/ data itrans/1, 2, 3, 7*0/ C C Predicted instrumental throughputs, taken from C the August 1995 Observers' Guide unless otherwise stated. C data effmirror/0.7/ C data effmirror/0.65/ ! 1st estim. (ARIST) C data effmirror/0.85/ data ((effinstr(i,j),i=1,10),j=1,10)/ C U B V R I Z J H K C IDS 235 throughput = 0.45 (camera) * 0.9 (collimator): C 1 0.41,0.41,0.41,0.41,0.41,0.41,0.0,0.0,0.0,0.0, ! INT IDS 235 camera C No optics: 1 1.00,1.00,1.00,1.00,1.00,1.00,0.0,0.0,0.0,0.0, ! JKT imaging 1 0.00,0.06,0.06,0.06,0.00,0.00,0.0,0.0,0.0,0.0, ! Human eye 1 80*0.0/ data (((effgrat(i,j,k),i=1,10),j=1,20),k=1,10)/ C 1 0.4,0.55,0.65,0.55,0.40,0.0,0.0,0.0,0.0,0.0, !IDS 235 C 1 0.4,0.55,0.7,0.65,0.4,0.0,0.0,0.0,0.0,0.0, C 1 0.55,0.45,0.3,0.0,0.0,0.0,0.0,0.0,0.0,0.0, C 1 0.4,0.65,0.65,0.5,0.35,0.0,0.0,0.0,0.0,0.0, C 1 0.0,0.0,0.4,0.55,0.55,0.55,0.0,0.0,0.0,0.0, C 1 0.0,0.0,0.55,0.6,0.55,0.5,0.0,0.0,0.0,0.0, C 1 0.0,0.0,0.0,0.3,0.65,0.7,0.0,0.0,0.0,0.0, C 1 0.4,0.6,0.65,0.6,0.45,0.4,0.0,0.0,0.0,0.0, C 1 0.0,0.0,0.4,0.55,0.65,0.65,0.0,0.0,0.0,0.0, C 1 0.0,0.65,0.65,0.55,0.0,0.0,0.0,0.0,0.0,0.0, C 1 0.7,0.6,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, C 1 0.65,0.6,0.5,0.4,0.0,0.0,0.0,0.0,0.0,0.0, C 1 0.35,0.55,0.65,0.65,0.4,0.4,0.0,0.0,0.0,0.0, C 1 0.0,0.0,0.4,0.6,0.65,0.0,0.0,0.0,0.0,0.0, C 1 0.0 ,0.3 ,0.45,0.5,0.0,0.0,0.0,0.0,0.0,0.0, C 1 0.25,0.55,0.5 ,0.0,0.0,0.0,0.0,0.0,0.0,0.0,40*0.0, 1 2000*1.0/ data ((dispmm(i,j),i=1,20),j=1,10)/ C 1 271.3,138.5,104.5,104.5,104.4,69.8,70.2,66.5, ! IDS235 C 1 50.7, 46.4, 35.3, 35.3, 35.2,34.8,23.2,17.5,4*0.0, 1 20*0.0, ! JKT imaging 1 180*0.0/ data arcmm/11.7, 13.8, 20000., 7*0.0/ data telarea/3.03,9*0.0/ data itel/1,1,2,7*0/ ! 1=Aristarchos, 2=eye data imirror/2, 3, 0, 7*0/ ! 2 for casegrain, 3 for all other instr, 0 for eye data itype/1, 2, 1, 7*0/ ! 1 for imaging, 2 for spectrum, 3 for fibre spect data ngrating/0, 3, 0, 7*0/ data iaxis/1, 1, 0, 7*0/ data idefdet /1, 2, 3, 7*0/ C 1 2 3 4 5 6 7 8 9 0 data defslit/1.0,1.0,1.0,7*1.0/ data spread /10*0.0/ C C Detector quantum efficiences etc. are from ING web page or articles C in SPECTRUM. Readout noise and gain are for typical readout C speeds, usually QUICK. For WHIRCAM, readout noise is for ND-STARE. C data effdet / 1 0.63,0.69,0.69,0.69,0.50,0.30,0.00,0.00,0.00,0.00, ! 1 TEK (UVAR - AR) C 1 0.30,0.60,0.65,0.68,0.45,0.25,0.00,0.00,0.00,0.00, ! 1 TEK 1 0.60,0.80,0.80,0.74,0.50,0.15,0.00,0.00,0.00,0.00, ! 2 EEVbig 1 0.00,1.00,1.00,1.00,0.00,0.00,0.00,0.00,0.00,0.00, ! 3 Eye 1 70*0.0/ data read /4.7,0.,0.,7*0.0/ C data read /6.,4.,0.,7*0.0/ data adu /4.15,0.8,1.0,7*0.0/ data pixel/0.024,0.0,0.003,7*0.0/ C data pixel/0.024,0.0135,0.003,7*0.0/ data npixels/1024,2048,6000,7*0, 1 1024,4096,6000,7*0/ data irot/0,1,0,7*0/ C C Labels: C data cband 1/'U ','B ','V ','R ','I ','Z ','J ','H ','Ks','? '/ data cbandlower 1/'u ','b ','v ','r ','i ','z ','j ','h ','ks','? '/ data (cinstr(i),i=1,10)/ 1 'ARISTARCHOS imaging ', 1 'INT IDS 235 camera ', 1 'Human eye (approx.) ', 1 7*' '/ data ((cgrating(i,j),i=1,20),j=1,10)/ 1 20*' ', ! JKT imaging 1 'R150V ','R300V ','R400B ', ! IDS 235 1 'R400V ','R400R ','R600R ', 1 'R600IR','R632V ','R831R ','R900V ','R1200U', 1 'R1200B','R1200Y','R1200R','H1800V','H2400B', 1 4*' ', ! IDS 235 1 160*' '/ data cdet(1) /'TEK ING 1024*1024, 78.0-mic'/ data cdet(2) /'EEV ING 2048*4096, 13.5-mic'/ C C Number of instruments, detectors and bands on offer C data ninstr/1/ data ndet /1/ data nband / 9/ C C Set default values of parameters C data iinstr/1/ data iband/3/ data idet/1/ data igrat/1/ data qgrat/' R158R'/ data fwhm/1/ data objm/20/ data slit/1/ data time/100/ data ipx/0/ data airmass/1.0/ C C############################################################################## C C Calculate PHOTONS from FLUXJY (v11.0 onward), set SKYMAG0 = SKYMAG C do i=1,nband photons(i)=fluxjy(i)*10000/6.6/wave(i) skymag0(i)=skymag(i) C print *,'Band, photons ',i,photons(i) enddo C C IP = 6 - interactive (normal mode) C 7 - non-interactive, reads data from a file signal.dat C writes output to signal.out C (e.g. for web interface) C 8 - non-interactive, reads data from a file standard input C writes output to standard output C (e.g. for web interface) C Non-interactive mode can also be triggered by typing X at the prompt. C C print *,'IP (6, 7 or 8)? ' C read(*,*) ip ip=6 C C If interactive, type welcome lines C if(ip.eq.6) then iin = 1 iout = 6 call welcome(version) else if (ip.eq.7) then iin = 1 iout = 7 open(unit=iin,name='/tmp/signal.dat') open(unit=iout,name='/tmp/signal.out') else iin = 5 iout = 6 endif C C Main menu loop C do loop=1,100000 C C Get parameters (from user or from SIGNAL.INPUT) C call getpar(ip,iin,iout,ninstr,ngrating,ndet,nband,loop, 1 cinstr,cgrating,cdet,cband,cbandlower, 1 effgrat, effdet, 1 wave,itel,idefdet,defslit,itrans, 1 iinstr,igrat, idet,iband, 1 objm,ipx,time,fwhm,slit,airmass,telarea,skymag,skymag0, 1 extin,band, 1 sol,eclat,glat,air,im,alpha,r,zm,ak, 1 mode) C C Set IPX = 0 for IFU/fibre spectroscopy C if(itype(iinstr).eq.3) ipx=0 C C Program flow now determined by MODE C C print *,'MODE = ',mode if(mode.eq.1) then call skybrcalc(ip,iout,skymag0, 1 sol,eclat,glat,air,im,alpha,r,zm,ak) goto 500 else if(mode.eq.2) then call help(ip,iout) goto 500 else if(mode.eq.3) then call helpsky(ip,iout) goto 500 else if(mode.eq.9) then goto 500 endif C C ELSE IF MODE=0, continue with S:N calculation C C C C If IP=7, write html tags C if(ip.eq.7) write(iout,492) 492 format('',/,'',/, 1'',/,'',/, 1'',/, 1'
') C C Convert input mag to required quantity C if(objm.gt.100) then objm=objm-100 objm=-1.e-6*10.**((23.9-objm)/2.5) write(iout,812) -objm 812 format('Brightness in Jy:',e10.3) endif if(objm.lt.0) then objm=-2.5*alog10(-objm/fluxjy(iband)) write(iout,811) objm 811 format('Apparent mag: ',f7.2) endif C C Reset QGRAT C if(igrat.ne.0) then qgrat=cgrating(igrat,iinstr)(1:10) else qgrat=' ' endif C C Calculate effective collecting area for TEL * INSTR * DET C geomarea=telarea(itel(iinstr)) area=geomarea 1 *10.**(-extin(iband)/2.5*airmass) 1 *effmirror**(float(imirror(iinstr))) 1 *effinstr(iband,iinstr)*effgrat(iband,igrat,iinstr) C print *,'area, geomarea = ',area,geomarea C print *,'extin, airmass, effmirror ',extin(iband),airmass, C 1 effmirror C print *,'imirror,effinstr,grat ',imirror(iinstr), C 1 effinstr(iband,iinstr),effgrat(iband,igrat,iinstr) throughput=area/geomarea area=area*effdet(iband,idet) if(area.eq.0) then write(iout,813) 813 format( 1 '*** WARNING - throughput = 0 Grating OK?', 1 ' Detector OK? Band OK?',/) endif C C For spectroscopic obs of point sources, calculate fraction C of light VIGN passing through slit, fibre or lenslet. C vign=1 if(itype(iinstr).ne.1.and.ipx.eq.0) 1 call slitvign(iinstr,itype,slit,fwhm,vign) C C Calculate scales at detector in arcsec/pixel and Angstrom/pixel C arcpix=arcmm(iinstr)*pixel(idet) disppix=dispmm(igrat,iinstr)*pixel(idet) C C Calculate number of photons from object and from sky/arcsec**2 C per Angstrom C sobj=photons(iband)*10000*area*time*10.**(-objm/2.5)*vign ssky= 1 photons(iband)*10000*area*time*10.**(-skymag(iband)/2.5) C C Correct for empirical/predicted throughput ratio C sobj=sobj*truethru(iband,iinstr) ssky=ssky*truethru(iband,iinstr) C C Display parameter settings to user C qinst=cinstr(iinstr) qband=cband(iband) qdet=cdet(idet)(1:6) qgrat2=qgrat sky=skymag(iband) deteff=effdet(iband,idet) pixsize=pixel(idet) photons2=photons(iband) bwidth=band(iband) readout=read(idet) slit2=slit true0=truethru(iband,iinstr) ext=extin(iband) if(ngrating(iinstr).gt.0) then qgrat2=qgrat else qgrat2=' ' endif if(itype(iinstr).eq.1) then slit2=0 else bwidth=0 endif call display(ip,iout,ans,iinstr,itype,ipx,ext, 1 iband,cband,wave, 1 qinst,qband,qgrat2,qdet,fwhm,objm,slit2,airmass, 1 sky,time,geomarea,throughput,deteff,pixsize,photons2, 1 bwidth,readout,true0,comment) C C Direct imaging calculation and output C if(itype(iinstr).eq.1) then write(iout,'('' '')') if(ipx.eq.0) then sobj2=sobj*band(iband) else sobj2=sobj*band(iband)*arcpix*arcpix endif ssky2=ssky*band(iband)*arcpix*arcpix call free(sobj2,buffa) call free(ssky2,buffb) call free(ssky,buffc) if(ipx.eq.0) then write(iout,821) buffa(1:10),buffc(1:10),buffb(1:10) 821 format('Detected photons from object =',a10,/, 1 'Detected photons/A/arcsec^2 from sky =',a10,/, 1 'Detected photons/pixel from sky =',a10) else if(ipx.eq.1) then write(iout,822) buffa(1:10),buffc(1:10),buffb(1:10) 822 format('Detected photons/pixel from object =',a10,/, 1 'Detected photons/A/arcsec^2 from sky =',a10,/, 1 'Detected photons/pixel from sky =',a10) endif else C C Spectroscopic calculation and output C call free(sobj,buffa) if(ipx.eq.0) then write(iout,823) buffa(1:10),vign 823 format('Detected photons/A from object =',a10, 1 ' (from in-slit fraction ',f4.2,')') else if(ipx.eq.1) then write(iout,824) buffa(1:10) 824 format( 1 'Detected photons per A from each sq arcsec of object =',a10) endif if(ipx.eq.0) then sobj2=sobj*disppix else sobj2=sobj*slit*arcpix*disppix endif if(sobj2.gt.0) then call free(sobj2,buffa) if(ipx.eq.0) then write(iout,825) buffa(1:10) 825 format( 1 'Detected photons/pixel-step-in-wavelength from object =', 1 a10) else if(ipx.eq.1) then write(iout,826) buffa(1:10) 826 format('Detected photons/pixel from object =',a10) endif endif ssky2=ssky*slit*arcpix*disppix if(itype(iinstr).eq.3) ssky2=ssky*slit*slit*3.14159/4.*disppix call free(ssky2,buffa) call free(ssky,buffc) if(ssky2.ne.0) then write(iout,734) buffc(1:10) 734 format('Detected photons/A/sq.arcsec from sky =',a10) if(itype(iinstr).eq.2) then write(iout,827) buffa(1:10) 827 format('Detected photons/pixel from sky =',a10) else if(itype(iinstr).eq.3) then write(iout,927) buffa(1:10) 927 format('Detected photons/pixel-step-in-wavelength', 1 ' from sky =',a10) endif endif endif C C Calculate signal-to-noise C C C Until version 11.0, used SEE = 3.14159/4.* for imaging C if(ipx.eq.0.and.arcpix.ne.0) then if(itype(iinstr).eq.1) see=3.14159*(fwhm/arcpix)**2 if(itype(iinstr).eq.2) see=2.*fwhm/arcpix else if(ipx.eq.1) then see=1. endif if(itype(iinstr).eq.3) see=spread(iinstr) if(itype(iinstr).le.2) then snr=sobj2/sqrt(sobj2+see*(ssky2+read(idet)*read(idet))) else if(itype(iinstr).eq.3) then snr=sobj2/sqrt(sobj2+ssky2+see*read(idet)*read(idet)) endif if(snr.ne.0) then if(arcpix.ne.0.or.itype(iinstr).eq.3) then call free(snr,buffa) if(ipx.eq.0) then if(itype(iinstr).eq.1) then write(iout,828) buffa(1:10) 828 format('Signal-to-noise =',a10) else if(itype(iinstr).eq.2) then write(iout,829) buffa(1:10) 829 format( 1 'Signal-to-noise/pixel of extracted spectrum =',a10) else if(itype(iinstr).eq.3) then write(iout,929) buffa(1:10) 929 format( 1 'Signal-to-noise/pixel of extracted spectrum', 1 ' (fibre or lenslet) =',a10) endif else if(ipx.eq.1) then write(iout,830) buffa(1:10) 830 format('Signal-to-noise/pixel = ',a10) endif endif endif C C Scales at the detector C For spectroscopy: C IAXIS(IINSTR) IROT(IDET) dispersion direction C 1 0 X C 2 0 Y C 1 1 Y C 2 1 X C write(iout,'('' '')') if(itype(iinstr).eq.1) then dx=npixels(idet,1)*arcpix dy=npixels(idet,2)*arcpix idx=1 idy=1 else if(iaxis(iinstr)+irot(idet).ne.2) then dx=npixels(idet,1)*disppix dy=npixels(idet,2)*arcpix idx=2 idy=1 else dx=npixels(idet,1)*arcpix dy=npixels(idet,2)*disppix idx=1 idy=2 endif endif C print *,'*** ',iinstr,iaxis(iinstr),idet,irot(idet),idx,idy if(dx.ne.0) write(iout,259) 1 dx/float(npixels(idet,1)),unit(idx), 1 npixels(idet,1),dx,unit(idx) 259 format('Scale in X direction: ',f7.2,1x,a1,'/pixel; ', 1 i4,' pixels = ',f9.2,1x,a1) if(dy.ne.0) write(iout,260) 1 dy/float(npixels(idet,2)),unit(idy), 1 npixels(idet,2),dy,unit(idy) 260 format('Scale in Y direction: ',f7.2,1x,a1,'/pixel; ', 1 i4,' pixels = ',f9.2,1x,a1) C C Photons/ADU conversion C write(iout,151) cdet(idet)(1:6),adu(idet) C ,int(read(idet)) 151 format(a6,' gain =',f4.1, 1 ' e-/ADU approx.') C ; readnoise =',i3,' electrons') C C If IP=7, write closing html tags C if(ip.eq.7) write(iout,493) 493 format('',/,'',/,'') enddo 500 continue end C****************************************************************************** C***** W E L C O M E ****************************************************** C****************************************************************************** C subroutine welcome(version) implicit none character*28 version print *,' ' print *,'SIGNAL '//version 1 //' crb@ing.iac.es' print *,' ' write(*,101) 101 format( 1'SIGNAL calculates photons/A detected from a source of given',/, 1'apparent magnitude (per arcsec**2 if extended), as:',/, 1' 10.**(-mag/2.5)',/, 1' * photons/sec/A/m**2 giving mag = 0 at top of atm.',/, 1' * transmission of atm. at given airmass',/, 1' * exposure time in sec',/, 1' * unobstructed area of main mirror in m**2',/, 1' * measured throughput of telescope/instrument',/, 1' * quantum efficiency of detector in the given band',/, 1'Accuracy +-20%. Signal-to-noise is per pixel step in',/, 1'wavelength for point sources,', 1' and per pixel for extended ones.',/, 1'For imaging, the counts obtained above are multiplied by the',/, 1'effective bandwidth of the filters in A, and signal-to-noise',/, 1'is calculated within a 2-FWHM-diameter aperture for a point',/, 1'source and per pixel for extended sources.') write(*,102) 102 format( 1'Signal-to-noise = Nobj/sqrt(Nobj+pixels*(Nsky+readnoise**2))',/, 1' ',/, 1'First time users, type H for more information.') end C****************************************************************************** C***** G E T P A R ******************************************************** C****************************************************************************** C subroutine getpar(ip,iin,iout,ninstr,ngrating,ndet,nband,loop, 1 cinstr,cgrating,cdet,cband,cbandlower, 1 effgrat, effdet, 1 wave,itel,idefdet,defslit,itrans, 1 iinstr,igrat, idet,iband, 1 objm,ipx,time,fwhm,slit,airmass,telarea,skymag,skymag0, 1 extin,band, 1 sol,eclat,glat,air,im,alpha,r,zm,ak, 1 mode) C C CMDOK string of acceptable user command characters C implicit none integer ninstr,ngrating(50),ndet,nband,itel(50),loop,iunit integer iinstr,igrat,idet,iband,ipx,mode,i,j,k,idefdet(50) integer itrans(50) real defslit(50) real effgrat(10,20,50),effdet(10,30),wave(10),objm,time real fwhm,slit,airmass,telarea(50),skymag(10),extin(10), 1 band(10),dummy character*1 ans character*2 cband(10),cbandlower(10),ans2 character*20 cgrating(20,50) character*25 cinstr(50) character*30 cinstrb(40) character*40 cdet(30) character*40 cmdok character*80 cline C C Variables used only in call to skybrcalc C integer ip,iin,iout,im real skymag0(10),sol,eclat,glat,air,alpha,r,zm,ak data cmdok/'IGDBMPTFSAKEWigdbmptfsakewQHqhUuXx '/ if(ip.ne.6) then C C ******************************** C C Get parameters from SIGNAL.INPUT C C ******************************** C mode=9 read(iin,*,end=500) mode if(mode.eq.0) then read(iin,*) iinstr read(iin,*) igrat if(igrat.eq.0) igrat=1 read(iin,*) idet read(iin,*) iband read(iin,*) objm read(iin,*) ipx read(iin,*) time read(iin,*) fwhm read(iin,*) slit read(iin,*) airmass read(iin,'(a80)') cline if(cline.eq.'D'.or.cline.eq.'d') then if(iband.le.6) skymag(iband)=skymag0(iband)-0.4 else if(cline.eq.'G'.or.cline.eq.'g') then if(iband.le.3) skymag(iband)=skymag0(iband)-2.15 if(iband.ge.4.and.iband.le.6) skymag(iband)=skymag0(iband)-1.3 else if(cline.eq.'B'.or.cline.eq.'b') then if(iband.le.3) skymag(iband)=skymag0(iband)-3.4 if(iband.ge.4.and.iband.le.6) skymag(iband)=skymag0(iband)-2.7 else read(cline,*) dummy if(dummy.ne.0) skymag(iband)=dummy endif read(iin,*) dummy if(dummy.ne.0) extin(iband)=dummy read(iin,*) dummy if(dummy.ne.0) band(iband)=dummy else if(mode.eq.1) then im=1 read(iin,*) sol read(iin,*) eclat read(iin,*) glat read(iin,*) air read(iin,*) alpha if(alpha.eq.-1) im=0 read(iin,*) zm read(iin,*) r read(iin,*) ak endif else C C ************************ C C Get parameters from user C C ************************ C iunit=5 print *,' ' 173 continue if(loop.gt.1) print *, 1 'Type I G D B M P T F S A U W K or E to change a parameter' print *, 1 'Type return to continue,', 1 ' Q to quit, H for more information' C print *,' ' read(iunit,'(a1)') ans C C Check command is one of set C do k=1,40 if(ans.eq.cmdok(k:k)) goto 174 enddo print *,'Invalid command: ',ans goto 173 174 continue C C Stop C if(ans.eq.'Q'.or.ans.eq.'q') mode=9 C C Type rudimentary help C if(ans.eq.'H'.or.ans.eq.'h') then call help(ip,iout) goto 27 endif C C Get band C if(ans.eq.'B'.or.ans.eq.'b') then write(*,193) 193 format('Approximate throughput for which of', 1 ' following bands?',/) do j=1,nband write(*,196) cband(j),int(wave(j)) 196 format(a2,' =',i6,' A for spectroscopy') enddo read(iunit,'(a2)') ans2 do j=1,nband if(ans2.eq.cband(j).or.ans2.eq.cbandlower(j)) iband=j enddo endif C C Get bandwidth C if(ans.eq.'W'.or.ans.eq.'w') then print *,'Current effective bandwidth for ',cband(iband) 1 ,' = ',band(iband) print *,'(= area under filter transmission curve)' print *,'New effective bandwidth for this band?' read(iunit,*) band(iband) endif C C Get apparent magnitude C if(ans.eq.'M'.or.ans.eq.'m') then write(6,104) cband(iband) 104 format(' Apparent magnitude in band ',a2,'?') print *, 1 '(number > 50 will be interpreted as 100 + Oke AB mag)' print *, 1 '(negative number will be interpreted as Jy', 1 ' (1E-26 W/Hz/m**2))' read(iunit,*) objm endif C C Get telescope/instrument combination C if(ans.eq.'I'.or.ans.eq.'i') then print *,'Telescope/instrument? (N/A = no longer offered)' do k=1,40 cinstrb(k)=' ' enddo do k=1,ninstr i=itrans(k) write(cinstrb(k),101) k,cinstr(i) 101 format(i3,2x,a25) enddo do k=1,20 write(6,'(a30,5x,a30)') cinstrb(k),cinstrb(k+20) enddo print *,' ' write(6,102) 102 format(' Choose a number ', 1 '(NB detector reverts to a default', 1 ' type on changing instrument)') read(iunit,*) k iinstr=itrans(k) idet=idefdet(iinstr) slit=defslit(iinstr) endif C C Get slit width C if(ans.eq.'S'.or.ans.eq.'s') then print *,'Slit width (or sqrt fibre/lenslet area) (arcsec)?' read(iunit,*) slit endif C C Get IPX (point or extended) C if(ans.eq.'P'.or.ans.eq.'p') then print *,'Mag is for point source (0) or', 1 ' per arcsec**2 for an extended source (1)?' read(iunit,*) ipx endif C C Get extinction C if(ans.eq.'E'.or.ans.eq.'e') then print *, 1 'New value of extinction per airmass for this band?' read(iunit,*) extin(iband) endif C C Get new value for telescope geometric area C if(ans.eq.'U'.or.ans.eq.'u') then print *,'Telescope geometric area = ', 1 telarea(itel(iinstr)) print *,'New value?' read(*,*) telarea(itel(iinstr)) endif C C Get airmass C if(ans.eq.'A'.or.ans.eq.'a') then print *,'Airmass?' read(iunit,*) airmass endif C C Get sky brightness C if(ans.eq.'K'.or.ans.eq.'k') then 729 print *,'Current value of sky brightness in ', 1 cband(iband),' = ',skymag(iband), 1 ' mags/arcsec**2' print *,'Enter new value for sky brightness' print *,'or enter 0 (zero) to run the sky-brightness', 1 ' calculator' print *,' (NB the calculator does not', 1 ' reset the current value)' print *,'or enter D, G or B to select', 1 ' typical sky brightness in', 1 ' dark, grey or bright of moon' print *,' (NB typical dark is not same as minimum)' print *,'or enter H for help' read(iunit,'(a80)') cline if(cline.eq.'H'.or.cline.eq.'h') then call helpsky(ip,iout) c goto 729 else if(cline.eq.'D'.or.cline.eq.'d') then if(iband.le.6) skymag(iband)=skymag0(iband)-0.4 else if(cline.eq.'G'.or.cline.eq.'g') then if(iband.le.3) skymag(iband)=skymag0(iband)-2.15 if(iband.ge.4.and.iband.le.6) skymag(iband)=skymag0(iband)-1.3 else if(cline.eq.'B'.or.cline.eq.'b') then if(iband.le.3) skymag(iband)=skymag0(iband)-3.4 if(iband.ge.4.and.iband.le.6) skymag(iband)=skymag0(iband)-2.7 else read(cline,*) dummy if(dummy.ne.0) then skymag(iband)=dummy else if(dummy.eq.0) then call skybrcalc(ip,iout,skymag0, 1 sol,eclat,glat,air,im,alpha,r,zm,ak) c goto 729 endif endif endif C C If NGRATING(IINSTR) > 0, get grating C if(ans.eq.'G'.or.ans.eq.'g') then print *, 1 'Grating? (NB for WYFFOS, the list below is', 1 ' 9 ISIS then 10 IDS gratings)' if(ngrating(iinstr).gt.0) then do i=1,ngrating(iinstr) write(6,105) i,cgrating(i,iinstr), 1 effgrat(iband,i,iinstr) 105 format(i2,2x,a20,' throughput in this band = ',f4.2) enddo write(6,172) 172 format('Choose a number (leftmost column above)') read(iunit,*) igrat endif endif C C Get detector C if(ans.eq.'D'.or.ans.eq.'d') then print *,'Detector?' print *, 1'(See ING web page for detailed detector info)' do i=1,ndet write(6,106) i,cdet(i),effdet(iband,i) 106 format(i3,2x,a40,' QE in this band = ',f4.2) enddo write(6,172) read(iunit,*) idet endif C C Get exposure time C if(ans.eq.'T'.or.ans.eq.'t') then print *,'Exposure time (sec)?' read(iunit,*) time endif C C Get FWHM C if(ans.eq.'F'.or.ans.eq.'f') then print *,'FWHM object * seeing?' read(iunit,*) fwhm endif 27 continue endif 500 continue end C****************************************************************************** C***** D I S P L A Y ****************************************************** C****************************************************************************** C subroutine display(ip,iout,ans,iinstr,itype, 1 ipx,ext, 1 iband,cband,wave, 1 qinst,qband,qgrat,qdet,fwhm,objm,slit,airmass, 1 sky,time,geomarea,throughput,deteff,pixsize,photons, 1 bwidth,readout,true0,comment) C C Displays parameter settings to user C implicit none integer i,iinstr,itype(50),ipx,iband,ip,iout real fwhm,objm,sky,time,geomarea,throughput,deteff,pixsize, 1 photons,bwidth,readout,slit,true0,ext,airmass,wave(10) character*6 qdet character*10 qgrat,qgratbuff character*37 text(20) character*25 qinst character*1 ans character*2 qband,cband(10) character*80 comment(50) do i=1,20 text(i)=' ' enddo qgratbuff=qgrat if(qgratbuff.eq.' ') qgratbuff=' none ' print *,' ' text(1)='G Grating '//qgratbuff text(2)='D Detector '//qdet C C Text strings which depend on whether imaging or spectroscopic C if(itype(iinstr).eq.1) then C C Imaging C text(3)='B Band '//qband write(text(4),'( 1 ''M Apparent magnitude '',f6.1)') objm write(text(16),'( 1 ''W Effective bandwidth (A) '',f9.2)') bwidth text(8)=' ' else C C Spectroscopic C if(itype(iinstr).eq.2) then write(text(8),'( 1 ''S Slit width (arcsec) '',f6.1)') slit else if(itype(iinstr).eq.3) then write(text(8),'( 1 ''S Fibre/lenslet diam (arcsec) '',f6.1)') slit endif if(ipx.eq.0) then write(text(3),'( 1 ''B Band '',a2, 1 '' = wavelength (A) '',i5)') 1 qband,int(wave(iband)) write(text(4),'( 1 ''M Apparent magnitude '',f6.1)') objm else write(text(4),'( 1 ''M Apparent magnitude/sq.arcsec'',f6.1)') objm endif text(15)=' ' endif C C Remaining text strings C write(text(5),'( 1 ''P Mag (0) or mag/sq.arcsec (1)'',i4)') ipx write(text(7),'( 1 ''F FWHM (object*seeing, arcsec)'',f6.1)') fwhm write(text(6),'( 1 ''T Integration time (sec) '',f9.1)') time write(text(9),'( 1 ''A Airmass '',f6.1)') airmass write(text(10),'( 1 ''U Tel. geometric area (sq.m) '',f6.2)') geomarea write(text(11),'( 1 '' Atm./tel./instr. throughput '',f6.2)') throughput write(text(12),'( 1 '' Detector efficiency '',f6.2)') deteff write(text(13),'( 1 '' Pixel size (microns) '',f6.1)') pixsize*1000 write(text(14),'( 1 '' Readout noise (electrons) '',f6.1)') readout write(text(15),'( 1 '' Photons/s/A/sq.cm mag=0 '',i5)') int(photons) write(text(18),'( 1 ''K Sky brightness mag/sq.arcsec '',f6.2)') sky write(text(19),'( 1 ''E Extinction per airmass (mag) '',f6.2)') ext write(text(17),'( 1 '' Empirical/theoretical '',f5.1)') true0 if(ip.eq.6) then write(iout,101) qinst,text(10),(text(i),text(i+10),i=1,9), 1 comment(iinstr) 101 format(' I Instrument = ',a25,1x,a37,/,8(1x,a37,4x,a37,/), 1 1x,a37,4x,a37,/,a80) else write(iout,102) qinst,text(10)(3:37), 1 (text(i)(3:37),text(i+10)(3:37),i=1,9), 1 comment(iinstr) 102 format('Instrument = ',a25,1x,a35,/,9(a35,4x,a35,/),a80) endif end C****************************************************************************** C***** S K Y B R C A L C ************************************************** C****************************************************************************** C subroutine skybrcalc(ip,iout,skymag, 1sol,eclat,glat,air,im,alpha,r,zm,ak) C C SKYBRCALC calculates the brightness of the sky in each band, C in S10 units and in mag/arcsec^2 as a function of C ecl, gal latitudes, zd, phase of solar cycle, and, if moon up: C phase of moon, zenith distance of moon, separation between C sky position and moon, extinction C If IP=6, the parameters are requested from user, otherwise C those supplied in call are used. C C Variables C --------- C AIR Airmass (derived from ZD) C AK V-band extinction in mag C ALPHA Phase of moon (0 = full, 90 = 7-day, 180 = new) C ECLAT Ec lat (deg) C GLAT Gal lat (deg) C IM = 1 if moonlight to be considered C IP Unit number to which to send output (6 = screen) C NBAND Number of wavelength bands C QAIR Airglow contribution in S10 C QSTAR Diffuse starlight contribution in S10 C QZOD Zodiacal light contibution in S10 C Q3 QAIR+QSTAT+QZOD in S10 C QSKY(10) Moonless sky brightness in S10 derived from SKY C QMOON(10) Contrib of moon in S10 C QALL(10) QSKY+QMOON C R Ang sep of patch of sky and moon (deg) C SKY Moonless sky br in mag/arcsec^2 derived from Q3, C correcting for colour difference (band - V) C SKYALL Total sky br in mag/arcsec^2 C SOL Phase of solar sycle (0 = min, 1 = max) converted t C 0.8 + 1.2*SOL C ZM Zenith distance of moon (deg) C implicit none integer i,im,nband,j,ip,iout,loop real skymag(10),sol,eclat,glat,air,alpha,r,zm,ak real qair,qzod,qstar,q3 real qsky(10),sky(10),qmoon(10),qall(10),skyall(10) nband=5 loop=1 if(ip.eq.6) loop=10000 do i=1,loop if(ip.eq.6) then print *,'--------------------------------------------------' print *,'Running sky-brightness calculator' print *, 1 'Solar activity (0 = minimum, 1 = maximum, 9 to finish', 1 ' with calculator)?' print *, 1 '(Dates of last three solar minima: 1976.8, 1986.8, 1996.5)' read(*,*) sol if(sol.eq.9) goto 22 print *, 1 'Ecliptic latitude in deg (longitude', 1 ' assumed to be > 90 deg from sun)?' read(*,*) eclat print *,'Galactic latitude in deg?' read(*,*) glat print *,'Zenith distance (deg)?' read(*,*) air air=1./sqrt(1-0.96*(sin(air/57.3))**2.) print *,'Phase angle of moon (-1 for moonless sky)?' print *,'(e.g. 0 = full, 90 = 7-day old moon, 180 = new moon)' read(*,*) alpha im=0 if(alpha.ne.-1) then im=1 print *,'Zenith distance of moon (deg)?' read(*,*) zm print *, 1 'Angular separation of moon and target patch of sky (deg)?' read(*,*) r print *,'V-band extinction (mag)?' read(*,*) ak endif endif ! if IP=6 C C Convert SOL to MJy C sol=0.8+1.2*sol call skybrcalc2(im,nband,skymag,sol,eclat,glat,air, 1 alpha,r,zm,ak, 1 qair,qzod,qstar,q3, 1 qsky,sky,qmoon,qall,skyall) C C If IP=7|8, write html header C if(ip.ne.6) write(iout,494) 494 format('',/,'',/, 1'',/,'',/, 1'',/, 1'Predicted brightness of La Palma sky', 1'',/,'') write(iout,116) (sol-0.8)/1.2,eclat,glat,acos(1/air)*57.3,air 116 format( 1 'Solar activity =',f4.1,', ecliptic and galactic', 1 ' latitudes =',f5.1,',',f5.1,' deg,',/, 1 ' ZD = ',f4.1,' deg, airmass = ',f5.3) if(im.eq.1) write(iout,117) alpha,zm,r,ak 117 format('Lunar phase = ',f5.1,' deg, lunar zenith distance =', 1 f5.1, ' deg,',/, 1 ' distance from moon = ',f5.1,' deg, extinction =', 1 f5.2,' mag') write(iout,101) 1 q3,qair,qzod,qstar, 1 skymag(1)-skymag(2),skymag(2)-skymag(3),skymag(3)-skymag(4), 1 skymag(4)-skymag(5), 1 (qsky(j),j=1,nband),(sky(j),j=1,nband) 101 format( 1 'Moonless V sky brightness, in S10 units ', 1 '(equivalent 10th-mag stars/deg^2):',/, 1 f6.1,' ( =',f6.1,' airglow +',f6.1,' zodiacal light +',f6.1, 1 ' starlight)',/, 1 'The (moonless) sky brightness in other bands is obtained by', 1 ' assuming ',/, 1 'colours U-B, B-V, V-R, R-I = ',4f6.2,':',//, 1 ' Band ', 1 'U B V R I',/, 1 ' Total (no moon)',5f9.1,' S10 units',/, 1 ' Total (no moon) ',5f9.2,' mag/arcsec^2') if(im.eq.1) 1 write(iout,102) 1 (qmoon(j),j=1,nband),(qall(j),j=1,nband), 1 (skyall(j),j=1,nband),(sky(j)-skyall(j),j=1,nband) 102 format(/, 1 ' Moonlight ',5f9.1,' S10 units',/, 1 ' Total ',5f9.1,' S10 units',/, 1 ' Total ',5f9.2,' mag/arcsec^2',/, 1 ' (Change with moon',5f9.2,' mag)',/) write(iout,109) 109 format('Accuracy:', 1 ' +- 0.1 mag without moon (in B,V),', 1 ' +- 1 mag with moon',/) C C If IP=7|8, write closing html tags C if(ip.ne.6) write(iout,495) 495 format('',/,'',/,'') enddo 22 continue end C****************************************************************************** C***** S K Y B R C A L C 2 ************************************************ C****************************************************************************** C subroutine skybrcalc2(im,nband,skymag,sol,eclat,glat,air, 1 alpha,r,zm,ak, 1 qair,qzod,qstar,q3, 1 qsky,sky,qmoon,qall,skyall) C C Calculates sky brightness for given conditions, using LP tech note C 115 (Benn & Ellison) for dark-sky brightness, and Krisciunas & Schaefer C 1991 PASP 103 1033 p.1038 for moonlight contribution. C All working in S10 units until final conversion to mag/arcsec**2 C C Variables not used in SKYBRCALC: C S, FR, XZ, XZM, BNL - terms at start Sect 5 of Krisciunas C & Schaefer for calculating moonlight contr. C BNL is moonlight sky brightness in V in nanoLamberts C FUDGE(I) normalise the Krisciunas prediction to values measured C in 6 bands at the JKT 5/7/98 C BS10 = BNL * 3.8, the V moonlight sky br in S10 units C implicit none integer i,im,nband real sol,eclat,glat,air,qair,qzod,qstar,q3 real qsky(10),sky(10),qmoon(10),qall(10),skyall(10) real skymag(10),fudge(10) real alpha,s,fr,r,xz,zm,xzm,bnl,ak,bs10 data fudge/2.0,1.6,2.4,3.0,5.9,5*0.0/ C C Calculate V sky brightness without moon C qair=(145+130*(sol-0.8)/1.2)*air if(eclat.lt.60) then qzod=140-90.*sin(eclat/57.3) else qzod=60 endif qstar=100.*exp(-abs(glat)/10.) q3=qair+qzod+qstar C C Assume sky colour implied by values in SKYMAG, and calculate C dark-sky brightness in the other bands C do i=1,nband sky(i)=27.78-2.5*log10(q3)+skymag(i)-skymag(3) qsky(i)=10.**((27.78-sky(i))/2.5) enddo C C Calculate contribution of moon in each band C if(im.eq.1) then s=10.**(-0.4*(3.84+0.026*alpha+4.e-9*alpha**4.)) fr=10.**5.36*(1.06+(cos(r/57.3))**2)+10.**(6.15-r/40.) xz =air xzm=1./sqrt(1-0.96*(sin(zm/57.3))**2.) bnl=s*fr*10.**(-0.4*ak*xzm)*(1-10.**(-0.4*ak*xz)) bs10=bnl*3.8 C print *,'s,fr,xz,xzm,bnl,bs10 ',s,fr,xz,xzm,bnl,bs10 do i=1,nband qmoon(i)=fudge(i)*bs10 qall(i)=qsky(i)+qmoon(i) skyall(i)=27.78-2.5*log10(qall(i)) enddo endif end C****************************************************************************** C***** H E L P ************************************************************ C****************************************************************************** C subroutine help(ip,iout) C C C implicit none integer ip,iout write(iout,'('' '')') write(iout,101) 101 format( 1'Full details of the program''s operation are given in the',/, 1'comments preceding the code.',/, 1'More up-to-date help is given on the SIGNAL web page:',/, 1' http://www.ing.iac.es/~crb/signal.html',/, 1' ',/, 1'For imaging, the bandwidths used correspond to the',/, 1'50-mm or liq-CuSO4 U filters, the Harris B,V,R,I set,',/, 1'the RGO glass Z filter and the WHIRCAM J, H and Kshort',/, 1'filters. Note that the 125-mm glass U filter',/, 1'has a factor 3 lower throughput than the 50-mm U.',/, 1'Note also that the throughputs of the Harris B filters are',/, 1'up to 25% less than that of the KPNO B filters.',/, 1'To specify narrow-band filters or adjust the effective',/, 1'bandwidth, use command B.',/, 1' ',/, 1'Object brightnesses can also be given in Jy or Oke AB mag',/, 1'(use the M command).') if(ip.eq.6) call pause write(iout,102) 102 format( 1' ',/, 1'The program does not take into account:',/, 1'- Vignetting at large field radius',/, 1'- Vignetting by the spectrograph slit (the specified width',/, 1' of the slit is used only in calculating the contribution',/, 1' from the sky).',/, 1'- Loss of light to colour and ND filters or to polarisation',/, 1' optics in spectrograph',/, 1'- Loss of grating efficiency at large angles',/, 1' and loss of light to neighbouring orders',/, 1'- Brightening of the sky due to twilight, moonlight,',/, 1' solar activity, low ecliptic or galactic latitude,',/, 1' or light pollution (enhanced in dusty conditions)',/, 1' - see option K in the menu',/, 1'- For FOS, U,B are for second order, V,R,I are first order.') if(ip.eq.6) call pause write(iout,103) 103 format( 1' ',/, 1'The signal-to-noise calculation:',/, 1' Nobj = photons/A (per arcsec**2 if extended) from object',/, 1' Nsky = photons/A/arcsec**2 from sky',/, 1' BAND = equivalent width of filter in A (integral T(l)dl',/, 1' where T(l) is transmission, l is wavelength)',/, 1' P = number of pixels over which integration carried out',/, 1' READ = CCD readout noise (e-)',/, 1' FWHM = object fwhm (intrinsic and due to seeing) in arcsec',/, 1' ',/, 1' Imaging, point or extended sources:',/, 1' N = Nobj * BAND for point sources',/, 1' N = Nobj * BAND*(arcsec/pixel)**2 for extended sources',/, 1' S = Nsky * BAND*(arcsec/pixel)**2',/, 1' P = pi*(FWHM/(arcsec/pixel))**2 for point source',/, 1' P = 1 for an extended source') if(ip.eq.6) call pause write(iout,104) 104 format( 1' ',/, 1' Spectroscopy, point source:',/, 1' N = Nobj * A/pixel',/, 1' S = Nsky * slit-width(A) * arcsec/pixel * A/pixel',/, 1' P = 2 * FWHM (arcsec) / arcsec/pixel',/, 1' ',/, 1' Spectroscopy, extended source:',/, 1' N = Nobj * slit-width(A) * arcsec/pixel * A/pixel',/, 1' S = Nsky * slit-width(A) * arcsec/pixel * A/pixel',/, 1' P = 1',/, 1' ',/, 1' Signal-to-noise = N/sqrt(N+P*(S+READ**2))',/, 1' For point sources, the sky counts are those from a',/, 1' 2-FWHM-diameter circular aperture for imaging, and from',/, 1' a 2-FWHM * slit-width rectangle for spectroscopy.') write(iout,105) 105 format( 1' For extended sources, the calculations are per pixel.') end C****************************************************************************** C***** H E L P S K Y ****************************************************** C****************************************************************************** C subroutine helpsky(ip,iout) implicit none integer ip,iout character*1 pans character*70 prompt pans=' ' write(iout,101) 101 format( 1'SIGNAL''s default sky-brightness settings are median for',/, 1'dark-of-moon, solar minimum, at high galactic and ecliptic',/, 1'latitude and in the absence of twilight and moonlight.',/, 1'The BVR values are accurate to +- 0.1 mag and are similar to',/, 1'those measured at other dark sites (Chile, Hawaii etc.).',/, 1'U and I are accurate to 0.5 mag. R and I sky brightnesses',/, 1'vary randomly by several tenths of a mag with variations in',/, 1'the OH airglow. The V and R sky brightness include a',/, 1'contribution of about 0.1 mag due to NaD light pollution.',/, 1'Light pollution is negligible in other bands, and it is',/, 1'decreasing.',/, 1'The sky is markedly brighter (several tenths of a mag) under',/, 1'very dusty conditions (> 0.3 mag extinction).',/, 1'The sky is 0.4 mag brighter at solar maximum. Recent minima') prompt= 1'[Return to continue, Q to end help]' if(ip.eq.6) call pause2(pans,prompt) if(pans.eq.'Q'.or.pans.eq.'q') goto 22 write(iout,102) 102 format( 1'were in 1986.5, 1996.5. Variation is approximately',/, 1'sinusoidal with time.',/, 1'The sky is 0.4 mag brighter on the ecliptic than at the',/, 1'poles, varying as sine(b) approx.',/, 1'The airglow contribution (typically about 70% of the total',/, 1'in V) brightens approx as airmass. The sky is 0.3 mag',/, 1'brighter at airmass 1.5.',/, 1'Stars fainter than mag=20 contribute', 1' negligibly to the total.',/, 1'Starlight scattered by interstellar dust contributes about',/, 1'5% of the total, rising to about 30% on the galactic plane.',/, 1'The extragalactic contribution is negligible (< 1%).',/, 1'The brightness of the sky does not vary with time after',/, 1'astronomical twilight.',/, 1'For further details of the calculations of moonless sky',/, 1'brightness, see La Palma technical note 115',/, 1'(Benn & Ellison 1998).') if(ip.eq.6) call pause2(pans,prompt) if(pans.eq.'Q'.or.pans.eq.'q') goto 22 write(iout,103) 103 format( 1' ',/, 1'SKY BRIGHTNESS WITH MOON UP:',/, 1'When the moon is 60 deg from zenith,', 1' with extinction 0.15 mag,',/, 1'the zenith sky will brighten roughly by M as tabulated below:',/, 1' ',/, 1' New Crescent Quarter Gibbous Full',/, 1' Phase angle (deg) 180 135 90 45 0 ',/, 1' Approx day: 1 4 8 12 15',/, 1' D, G or B: D G G B B',/, 1' Illum. frac. % 0 25 50 75 100',/, 1' M (U, B, V) 0 0.5 2.0 3.1 4.3',/, 1' M (R) 0 0.3 1.3 2.4 3.5',/, 1' M (I) 0 0.2 1.1 2.2 3.3',/, 1' ') if(ip.eq.6) call pause2(pans,prompt) if(pans.eq.'Q'.or.pans.eq.'q') goto 22 write(iout,104) 104 format( 1'Sky brightness for other values of lunar phase, lunar zenith',/, 1'angle, sky position and extinction, can be estimated with',/, 1'SIGNAL''s sky-brightness calculator.',/, 1'Note that the quarter moon (i.e. half disc illuminated)',/, 1'is a factor of 10 (not 2) fainter than full, due to the',/, 1'opposition effect (also responsible for gegenschein on the',/, 1'ecliptic and dry heiligenschein on earth).',/, 1'The contribution of moonlight in V has been calculated',/, 1'according to the scattering formula of Krisciunas & Schaefer',/, 1'(1991, PASP, 103, 1033), normalised (multiplied by a factor',/, 1'of 2.4) to agree with measurements of sky brightness made at',/, 1'the JKT on a dust-free night in 7/98.',/, 1'The moonlight contribution in the other bands is calculated',/, 1'according to the U-B, B-V, V-R, R-I colours of moonlight',/, 1'measured on the same night in 7/98.') if(ip.eq.6) call pause2(pans,prompt) if(pans.eq.'Q'.or.pans.eq.'q') goto 22 write(iout,105) 105 format( 1'These values agree +-40% with measurements made by DHPJ in',/, 1'9/89, but the contribution of moonlight probably depends',/, 1'strongly on local conditions (e.g. dust, telescope baffling),',/, 1'and with current data, the contribution by moonlight on',/, 1'La Palma can probably only be predicted within a factor ~2.',/, 1' ',/, 1'Ian Steele (LJMU) has found that the background brightness',/, 1'at the JKT rises dramatically (factor >~5 brighter than',/, 1'the above numbers) if moonlight falls on the telescope',/, 1'structure (scattering within the telescope).') 22 continue end C****************************************************************************** C***** P A U S E ********************************************************** C****************************************************************************** C subroutine pause implicit none character*1 ans print *,'[Press return to continue]' read(*,'(a1)') ans end subroutine pause2(pans,prompt) implicit none character*1 pans character*70 prompt print *,prompt read(*,'(a1)') pans end subroutine free(x,buff15) C C Writes real X into buffer BUFF15 free-format, as in PRINT *,X C (unix fortran print * not convenient because renders numbers C outside 0.1049 - 1012863 (why?) in exponential format, and may C not yield same formats as other fortran) C implicit none real x character*15 buff15 C print *,'okfree 1',x,buff15 if(x.lt.1000000.and.x.gt.0.1) then write(buff15,'(f10.2,5x)') x else write(buff15,'(e10.3,5x)') x endif C print *,'okfree 2',x,buff15 end subroutine slitvign(iinstr,itype,slit,fwhm,vign) C C Calculates VIGN, the fraction of light entering slit for given C pointlike object FWHM. C implicit none integer iinstr,itype(50) real slit,fwhm,vign,sigma,r,rat vign=1 C C Slit spectroscopy C Formulae below are approx (better than 5%) C fit to dependence of light fraction in slit vs slit/fwhm C given by light_in_slit.f numerical simulations C if(itype(iinstr).eq.2) then if(fwhm.eq.0) then vign=1.0 else rat=slit/fwhm if(rat.lt.0.76) vign=0.868*rat if(rat.ge.0.76.and.rat.lt.1.40) vign=0.37+0.393*rat if(rat.ge.1.40.and.rat.lt.2.30) vign=1.00-0.089*(2.3-rat) if(rat.ge.2.3) vign=1.0 endif endif C C Fibre/lenslet spectroscopy C r=slit/2. sigma=fwhm/2.35 if(itype(iinstr).eq.3) then vign=1.0-exp(-r*r/2./sigma/sigma) endif end