C @(#)airmass.for 17.1.1.1 (ES0-DMD) 01/25/02 17:40:55 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM AIRMAS C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.NAME C AIRMASS.FOR VERSION 3.0 C C.AUTHOR C M. ROSA (STECF) VERSION 1.0 29-OCT-1985 C C.MODIFICATIONS C M. ROSA (STECF) VERSION 2.0 29-MAR-1986 C Add JD and Sid. Time part, correct bug in airmass part C Input is in DEGREES for RA, DEC, Long., Lat. C Add more output C M. ROSA (STECF) VERSION 3.0 30-MAY-1986 C Verify for inclusion in MIDAS release C D. Baade (ST-ECF) VERSION 4.0 07-FEB-1987 C a) chop off calculation of ST and JD C b) improved formula for calculation of airmass as a function of C zenith distance C D. Baade (ESO) VERSION 4.1 21-DEC-1988 C Use more accurate weights for calculation of mean airmass C in the case of long exposures (formula given by P. Stetson, C 1988 September, Dominion Astrophysical Observatory Preprint) C K. Banse 890213 C move to new ST interfaces and FORTRAN 77 C C.PURPOSE a) Calculate airmass through sec z C C Correct expansion C C Airm = sec(z) - 0.0018167(sec(z)-1) C - 0.002875 (sec(z)-1)**2 C - 0.0008083(sec(z)-1)**3 C C.INPUT C hour angle (real degrees) INPUTR(1:1) C declination (degree,min,sec) INPUTR(2:4) C exposure time (sec) INPUTR(5:5) C latitude (degree,min,sec) INPUTR(6:8) C C.OUTPUT C OUTPUTR(1) Airmass C C--------------------------------------------------------------- C IMPLICIT NONE C REAL RB(8),AIR,DELT,WEIGHT(3) C DOUBLE PRECISION HA,DD,SECZ,LAT DOUBLE PRECISION FACTO,PI,RLAT,RDD,RHA C INTEGER I,IACTV,ISTAT,KUN(1),KNUL C CHARACTER*2 SIGNS C DATA PI /3.1415926535D0/ DATA WEIGHT/1.,4.,1./ C CALL STSPRO('AIRMAS') C CALL STKRDC('SIGNS',1,1,2,IACTV,SIGNS,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,8,IACTV,RB,KUN,KNUL,ISTAT) C FACTO = PI / 180.0D0 C HA = RB(1) ! hour angle in degrees 10 IF (HA.LT.-180.0D0) THEN HA = HA+360. GOTO 10 ENDIF C 11 IF (HA.GT.180.0D0) THEN HA = HA-360. GOTO 11 ENDIF C C declination in real degrees DD = ABS(RB(2))+(ABS(RB(3))+ABS(RB(4))/60.)/60. IF (SIGNS(2:2).EQ.'-') DD = -DD !negative declination C C latitude in real degrees LAT = ABS(RB(6))+(ABS(RB(7))+ABS(RB(8))/60.)/60. IF (SIGNS(1:1).EQ.'-') LAT = -LAT !negative latitude RLAT = LAT * FACTO !move from degrees to radians RDD = DD * FACTO C C first case: exposure time infinitly short C IF (RB(5).LE.0.0) THEN RHA = HA * FACTO SECZ = DSIN(RLAT)*DSIN(RDD)+ + DCOS(RLAT)*DCOS(RDD)*DCOS(RHA) SECZ = (1.0D0/SECZ) ! zenith distance C C calulate airmass C (formula by Young and Irvine, Astron. J. 72, 945 (1967) used below C includes refraction correction but is in principle valid only C for observations made at sea level) C AIR = SNGL(SECZ*(1.-0.0012*(SECZ*SECZ-1.))) C C Second case: exposure time > 0.0. Then, the airmass calculated is the C weighted average of airmasses at t0, 0.5(t1-t0) and t1 with weighting C factors w0 = 1/6, w1/2 = 4/6, and w1 = 1/6: C C W A R N I N G: What follows is a linear interpolation of a C hyperexponentially varying function ! C ELSE DELT = RB(5)/480.0D0 ! 0.5 * (sec)*15/3600 HA = HA-DELT AIR = 0.0 C DO 1000, I=1,3 HA = HA+DELT RHA = HA * FACTO SECZ = DSIN(RLAT)*DSIN(RDD)+ + DCOS(RLAT)*DCOS(RDD)*DCOS(RHA) SECZ = (1.0D0/SECZ) ! zenith distance C C calulate airmass C (below formula by Young and Irvine, Astron. J. 72, 945 (1967) C includes refraction correction but is in principle valid only C for observations made at sea level) C AIR = AIR + WEIGHT(I)*SNGL(SECZ*(1.-0.0012*(SECZ*SECZ-1.))) 1000 CONTINUE C AIR = AIR/6.0 ENDIF C C write airmass to keyword OUTPUTR CALL STKWRR('OUTPUTR',AIR,1,1,KUN,ISTAT) CALL STSEPI C END