C @(#)dvecs.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:33
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
      SUBROUTINE DOPFC(OPCODE,A,SCALR,C,NDIM)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine DOPFC            version 2.00      880414
C             DOPFF                    2.00      880414
C  K. Banse                        ESO - Garching
C  source code in file dvecs.for
C  2.10      890918
C
C.KEYWORDS
C  double precision, bulk data frames, arithmetic operations
C
C.PURPOSE
C  perform arithmetic on arrays
C
C.ALGORITHM
C  do it according to opcode ( `, +, -, *, / or =)
C  ` means ** (to the power of)
C  special condition handler is used for division by "zero"
C
C.INPUT/OUTPUT
C  call as OPFC(OPCODE,A,SCALR,C,NDIM)
C          OPFF(OPCODE,A,B,C,NDIM)
C
C  input par:
C  OPCODE:      char.exp.      operation code in postfix notation, e.g. "FF+"
C  A:           R*8 array      1. input array
C  B:           R*8 array      2. input array
C  SCALR:       R*8            scalar value
C  CONSTA:      R*8            1. input scalar
C  CONSTB:      R*8            2. input scalar
C  NDIM:        I*4            size of arrays involved
C
C  output par:
C  C:           R*8 array      result array
C  CONSTC:      R*8            result scalar
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NDIM,N,NCOUNT,TSTFLG
C
      DOUBLE PRECISION    A(NDIM),C(NDIM),SCALR,SCALRR
      DOUBLE PRECISION    USRNUL,DEPSP,DEPSN
      DOUBLE PRECISION    ONEG,ONEL
C
      CHARACTER*(*)      OPCODE
C
      COMMON      /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN
C
      DATA        ONEG  /1.000000001/
      DATA        ONEL  /0.999999999/
C 
C  subtract 2. operand from 1. operand
C 
      IF (OPCODE(3:3).EQ.'-') THEN
         IF (OPCODE(1:1).EQ.'C') THEN
            IF ((SCALR.LT.DEPSN) .OR. (SCALR.GT.DEPSP)) THEN
               DO 200, N=1,NDIM
                  C(N) = SCALR - A(N)
200            CONTINUE
            ELSE                                  !0. - C(n) 
               DO 220, N=1,NDIM
                  C(N) = -A(N) 
220            CONTINUE
            ENDIF
         ELSE
            IF ((SCALR.LT.DEPSN) .OR. (SCALR.GT.DEPSP)) THEN
               DO 300, N=1,NDIM
                  C(N) = A(N) - SCALR
300            CONTINUE
            ELSE                                  !C(n) - 0. 
               DO 330, N=1,NDIM
                  C(N) = A(N) 
330            CONTINUE
            ENDIF
         ENDIF
C
C  add operands
C 
      ELSE IF (OPCODE(3:3).EQ.'+') THEN
         IF ((SCALR.LT.DEPSN) .OR. (SCALR.GT.DEPSP)) THEN
            DO 1000, N=1,NDIM
               C(N) = A(N) + SCALR
1000        CONTINUE
         ELSE                                  !C(n) - 0. 
            DO 1010, N=1,NDIM
               C(N) = A(N) 
1010        CONTINUE
         ENDIF
C
C  divide 1. operand by 2. operand
C 
      ELSE IF (OPCODE(3:3).EQ.'/') THEN
         IF (OPCODE(1:1).EQ.'C') THEN
            IF (TSTFLG.EQ.1) THEN
               DO 1300, N=1,NDIM
                  IF ((A(N).LT.DEPSP).AND. (A(N).GT.DEPSN)) THEN
                     C(N) = USRNUL
                     NCOUNT = NCOUNT + 1
                  ELSE
                     C(N) = SCALR / A(N)
                  ENDIF
1300           CONTINUE
            ELSE
               DO 1350, N=1,NDIM
                  C(N) = SCALR / A(N)
1350           CONTINUE
            ENDIF
C
         ELSE                                        !divisor always tested 
            SCALRR = ABS(SCALR)
            IF ((SCALRR.LT.ONEG) .AND. (SCALRR.GT.ONEL)) THEN
               IF (SCALR.GT.0) THEN                  !C(n) / 1.
                  DO 1500, N=1,NDIM
                     C(N) = A(N) 
1500              CONTINUE
               ELSE                                  !C(n) / -1.
                  DO 1520, N=1,NDIM
                     C(N) = -A(N) 
1520              CONTINUE
               ENDIF
            ELSE IF (SCALRR.LT.DEPSP) THEN
               DO 1550, N=1,NDIM
                  C(N) = USRNUL
1550           CONTINUE
               NCOUNT = NCOUNT + NDIM
            ELSE
               SCALRR = 1.0/SCALR
               DO 1580, N=1,NDIM
                  C(N) = A(N) * SCALRR
1580           CONTINUE
            ENDIF
         ENDIF
C
C  multiply operands
C 
      ELSE IF (OPCODE(3:3).EQ.'*') THEN
         SCALRR = ABS(SCALR)
         IF ((SCALRR.LT.ONEG) .AND. (SCALRR.GT.ONEL)) THEN
            IF (SCALR.GT.0) THEN                  !C(n) * 1.
               DO 4000, N=1,NDIM
                  C(N) = A(N)
4000           CONTINUE
            ELSE                                  !C(n) * -1.
               DO 4020, N=1,NDIM
                  C(N) = -A(N)
4020           CONTINUE
            ENDIF
         ELSE IF (SCALRR.LT.DEPSP) THEN            !C(n) * 0.
            DO 4050, N=1,NDIM
               C(N) = 0.
4050        CONTINUE
         ELSE
            DO 4100, N=1,NDIM
               C(N) = A(N) * SCALR
4100        CONTINUE
         ENDIF
C
C  copy constant
C 
      ELSE IF (OPCODE(3:3).EQ.'=') THEN
         DO 4400, N=1,NDIM
            C(N) = SCALR
4400     CONTINUE
C
C  compute 1. operand ** 2. operand
C 
      ELSE IF (OPCODE(3:3).EQ.'`') THEN
         IF (OPCODE(1:1).EQ.'C') THEN
            DO 4800, N=1,NDIM
               C(N) = SCALR ** (A(N))
4800        CONTINUE
         ELSE
            IF ((SCALR.LT.DEPSN) .OR. (SCALR.GT.DEPSP)) THEN
               DO 5000, N=1,NDIM
                  C(N) = A(N) ** SCALR
5000           CONTINUE
            ELSE                                  !C(n) = x**0
               DO 5010, N=1,NDIM
                  C(N) = 1.
5010           CONTINUE
            ENDIF
         ENDIF
      ENDIF
C
      RETURN
      END

      SUBROUTINE DOPFF(OPCODE,A,B,C,NDIM)
C
      IMPLICIT NONE
C
      INTEGER      NDIM,N,NCOUNT,TSTFLG
C
      CHARACTER*(*)      OPCODE
C
      DOUBLE PRECISION   A(*),B(*),C(*)
      DOUBLE PRECISION   USRNUL,DEPSP,DEPSN
C
      COMMON      /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN
C
C  add frames
      IF (OPCODE(3:3).EQ.'+') THEN
         DO 200, N=1,NDIM
            C(N) = A(N) + B(N)
200      CONTINUE
C
C  subtract 2. frame from 1. frame
      ELSE IF (OPCODE(3:3).EQ.'-') THEN
         DO 400, N=1,NDIM
            C(N) = A(N) - B(N)
400      CONTINUE
C
C  multiply frames
      ELSE IF (OPCODE(3:3).EQ.'*') THEN
         DO 600, N=1,NDIM
            C(N) = A(N) * B(N)
600      CONTINUE
C
C  divide 1. frame by 2. frame
      ELSE IF (OPCODE(3:3).EQ.'/') THEN
         IF (TSTFLG.EQ.1) THEN
            DO 800, N=1,NDIM
               IF ((B(N).LT.DEPSP).AND. (B(N).GT.DEPSN)) THEN
                  C(N) = USRNUL
                  NCOUNT = NCOUNT + 1
               ELSE
                  C(N) = A(N) / B(N)      
               ENDIF
800         CONTINUE
         ELSE
            DO 850, N=1,NDIM
               C(N) = A(N) / B(N)      
850         CONTINUE
         ENDIF
C
C  copy 1. frame into result frame
      ELSE IF (OPCODE(3:3).EQ.'=') THEN
         DO 1000, N=1,NDIM
            C(N) = A(N)
1000     CONTINUE
C
C  compute 1. frame ** 2. frame
      ELSE IF (OPCODE(3:3).EQ.'`') THEN
         DO 1200, N=1,NDIM
            C(N) = A(N) ** B(N)
1200     CONTINUE
      ENDIF
C
      RETURN
      END

      SUBROUTINE DF1F(CFUNC,A,C,NDIM)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine DF1F            version 3.00      850822
C             DF1C            	3.00	850821
C             DF2FF            	3.00	850821
C             DF2FC            	3.00	850821
C             DF2CC            	3.00	850821
C  3.10	  910723
C 
C  K. Banse                        ESO - Garching
C
C.KEYWORDS
C  bulk data frames, functions
C
C.PURPOSE
C  apply usual FORTRAN functions to frames and constants
C
C.ALGORITHM
C  use special condition handler NULLER to trap invalid input to functions
C  all trigonometric functions work with degrees...
C
C.INPUT/OUTPUT
C  call as DF1F(CFUNC,A,C,NDIM)
C          DF1C(CFUNC,AVAL,CVAL)
C          DF2FF(CFUNC,A,B,C,NDIM)
C          DF2FC(CFUNC,A,SCALR,C,NDIM)
C          DF2CC(CFUNC,AVAL,CVAL)
C
C  input par:
C  CFUNC:      char.exp.      input function
C  NDIM:       I*4            size of arrays involved
C  A:          R*8 array      1. input array
C  B:          R*8 array      2. input array
C  AVAL:       R*8            input scalar value
C                  	input scalar values (2) for CCFUNC ...!
C
C  output par:
C  CVAL:       R*8            output scalar value
C  C:          R*8 array      result array
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER     N,NDIM,NCOUNT,TSTFLG
C
      DOUBLE PRECISION  A(NDIM),C(NDIM),FACT,FFACT
      DOUBLE PRECISION  USRNUL,DEPSP,DEPSN
C
      CHARACTER*5   CFUNC
C
      COMMON      /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN
C
      DATA  FACT  /0.017453292519444444/                  ! Pi / 180
      DATA  FFACT /57.2957795147199533/                   ! 180 / Pi 
C
C  branch according to function
C
C  SQRT
      IF (CFUNC.EQ.'SQRT ') THEN
         IF (TSTFLG.EQ.1) THEN
            DO 110, N=1,NDIM
               IF (A(N).LT.0.) THEN
                  C(N) = USRNUL
                  NCOUNT = NCOUNT + 1
               ELSE
                  C(N) = DSQRT(A(N))
               ENDIF
110         CONTINUE
         ELSE
            DO 120, N=1,NDIM
               C(N) = DSQRT(A(N))
120         CONTINUE
         ENDIF
C
C  LN
      ELSE IF (CFUNC.EQ.'LN   ') THEN
         IF (TSTFLG.EQ.1) THEN
            DO 210, N=1,NDIM
               IF (A(N).LE.0.) THEN
                  C(N) = USRNUL
                  NCOUNT = NCOUNT + 1
               ELSE
                  C(N) = DLOG(A(N))
               ENDIF
210         CONTINUE
         ELSE
            DO 220, N=1,NDIM
               C(N) = DLOG(A(N))
220         CONTINUE
         ENDIF
C
C  LOG10 or LOG
      ELSE IF (CFUNC(1:3).EQ.'LOG') THEN
         IF (TSTFLG.EQ.1) THEN
            DO 310, N=1,NDIM
               IF (A(N).LE.0.) THEN
                  C(N) = USRNUL
                  NCOUNT = NCOUNT + 1
               ELSE
                  C(N) = DLOG10(A(N))
               ENDIF
310         CONTINUE
         ELSE
            DO 320, N=1,NDIM
               C(N) = DLOG10(A(N))
320         CONTINUE
         ENDIF
C
C  EXP
      ELSE IF (CFUNC.EQ.'EXP  ') THEN
         DO 410, N=1,NDIM
            C(N) = DEXP(A(N))
410      CONTINUE
C
C  EXP10
      ELSE IF (CFUNC.EQ.'EXP10') THEN
         DO 510, N=1,NDIM
         C(N) = 10.D0**(A(N))
510   CONTINUE
C
C  SINE
      ELSE IF (CFUNC.EQ.'SIN  ') THEN
         DO 610, N=1,NDIM
         C(N) = DSIN ( A(N) * FACT )
610   CONTINUE
C
C  COSINE
      ELSE IF (CFUNC.EQ.'COS  ') THEN
         DO 710, N=1,NDIM
         C(N) = DCOS(A(N)*FACT)
710   CONTINUE
C
C  TANGENT
      ELSE IF (CFUNC.EQ.'TAN  ') THEN
         DO 810, N=1,NDIM
         C(N) = DTAN(A(N)*FACT)
810   CONTINUE
C
C  ARC SINE
      ELSE IF (CFUNC.EQ.'ASIN ') THEN
         DO 910, N=1,NDIM
            C(N) = DASIN(A(N)) * FFACT
910      CONTINUE
C
C  ARC COSINE
      ELSE IF (CFUNC.EQ.'ACOS ') THEN
         DO 1010, N=1,NDIM
            C(N) = DACOS(A(N)) * FFACT
1010     CONTINUE
C
C  ARC TANGENT
      ELSE IF (CFUNC.EQ.'ATAN ') THEN
         DO 1110, N=1,NDIM
            C(N) = DATAN(A(N)) * FFACT
1110     CONTINUE
C
C  NINT
      ELSE IF (CFUNC.EQ.'INT  ') THEN
         DO 1210, N=1,NDIM
            C(N) = DBLE(NINT(A(N)))
1210     CONTINUE
C
C  ABS
      ELSE IF (CFUNC.EQ.'ABS  ') THEN
         DO 1310, N=1,NDIM
            C(N) = DABS(A(N))
1310     CONTINUE
      ENDIF
C
      RETURN
      END

      SUBROUTINE DF2FF(CFUNC,A,B,C,NDIM)
C
      IMPLICIT NONE
C
      INTEGER      N,NDIM
C
      CHARACTER*(*)  CFUNC
C
      DOUBLE PRECISION   A(*),B(*),C(*)
      DOUBLE PRECISION   RA,RB,FACT,FFACT
C
      DATA  FACT  /0.017453292519444444/                  ! Pi / 180
      DATA  FFACT /57.2957795147199533/                   ! 180 / Pi 
C
C  branch according to function
C
C  function ATAN2
      IF (CFUNC.EQ.'ATAN2') THEN
         DO 1100, N=1,NDIM
            RA = A(N) * FACT
            RB = B(N) * FACT
            C(N) = DATAN2(RA,RB) * FFACT
1100     CONTINUE
C
C  minimum of 2 frames
      ELSE IF (CFUNC.EQ.'MIN  ') THEN
         DO 2100, N=1,NDIM
            C(N) = MIN(A(N),B(N))
2100     CONTINUE
C
C  maximum of 2 frames
      ELSE IF (CFUNC.EQ.'MAX  ') THEN
         DO 3100, N=1,NDIM
            C(N) = MAX(A(N),B(N))
3100     CONTINUE
C
C  remainder of 1.frame/2. frames
      ELSE IF (CFUNC.EQ.'MOD  ') THEN
         DO 4100, N=1,NDIM
            C(N) = DMOD(A(N),B(N))
4100     CONTINUE
      ENDIF
C
      RETURN
      END

      SUBROUTINE DF2FC(CFUNC,A,SCALR,C,NDIM)
C
      IMPLICIT NONE
C
      INTEGER    N,NDIM
C
      DOUBLE PRECISION  A(*),SCALR,C(*)
      DOUBLE PRECISION  RA,RB,FACT,FFACT
C
      CHARACTER*(*)   CFUNC
C
      DATA  FACT  /0.017453292519444444/                  ! Pi / 180
      DATA  FFACT /57.2957795147199533/                   ! 180 / Pi 
C
C  branch according to function
C
C  function ATAN2
      IF (CFUNC.EQ.'ATAN2') THEN
         RB = SCALR * FACT
         DO 1100, N=1,NDIM
            RA = A(N) * FACT
            C(N) = DATAN2(RA,RB) * FFACT
1100     CONTINUE
      RETURN
C
C  get minimum
      ELSE IF (CFUNC.EQ.'MIN  ') THEN
         DO 2100, N=1,NDIM
            C(N) = MIN(A(N),SCALR)
2100     CONTINUE
      RETURN
C
C  get maximum
      ELSE IF (CFUNC.EQ.'MAX  ') THEN
         DO 3100, N=1,NDIM
            C(N) = MAX(A(N),SCALR)
3100     CONTINUE
      RETURN
C
C  get remainder
      ELSE IF (CFUNC.EQ.'MOD  ') THEN
         DO 4100, N=1,NDIM
            C(N) = DMOD(A(N),SCALR)
4100     CONTINUE
      ENDIF
C
      RETURN
      END