C @(#)airmass.for	19.1 (ESO-DMD) 02/25/03 14:01:40
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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
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 020503	last modif	KB
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./
      DATA RB /8*0.0/,  SIGNS /'  '/
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
