C @(#)dcalc.for	19.2 (ESO-DMD) 05/20/03 09:41:54
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 DOPCC(OPCODE,CONSTA,CONSTB,CONSTC)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine OPCC           version 2.00      880202
C  K. Banse                  ESO - Garching
C  source code on file CALC.F
C  2.10      890918
C
C.KEYWORDS
C  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 OPCC(OPCODE,CONSTA,CONSTB,CONSTC)
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      NCOUNT
C
      DOUBLE PRECISION     CONSTA,CONSTB,CONSTC
      DOUBLE PRECISION     DEPSP,DEPSN
C
      REAL     USRNUL
C 
      CHARACTER*(*)      OPCODE
C
      COMMON      /NULCOM/ NCOUNT,USRNUL
C
      DATA DEPSP /10.D-33/
      DATA DEPSN /-10.D-33/
C
C  branch according to opcode
C
C  subtract 2. scalar from 1. scalar
      IF (OPCODE(3:3).EQ.'-') THEN
         CONSTC = CONSTA - CONSTB
C
C  add scalars
      ELSE IF (OPCODE(3:3).EQ.'+') THEN
         CONSTC = CONSTA + CONSTB
C
C  divide 1. scalar by 2. scalar
      ELSE IF (OPCODE(3:3).EQ.'/') THEN
         IF ((CONSTB.LT.DEPSP) .AND.
     +       (CONSTB.GT.DEPSN)) THEN		!always check divisor...
            CONSTC = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CONSTC = CONSTA / CONSTB
         ENDIF
C
C  multiply scalars
      ELSE IF (OPCODE(3:3).EQ.'*') THEN
         CONSTC = CONSTA * CONSTB
C
C  compute 1. scalar ** 2. scalar
      ELSE IF (OPCODE(3:3).EQ.'`') THEN
         CONSTC = CONSTA ** CONSTB
      ENDIF
C
      RETURN
      END

      SUBROUTINE DF1C(CFUNC,AVAL,CVAL)
C
      IMPLICIT NONE
C
      INTEGER    NCOUNT
C
      DOUBLE PRECISION       AVAL,CVAL
      DOUBLE PRECISION       TEMP,PI
C
      REAL     USRNUL
C
      CHARACTER*5  CFUNC
C
      COMMON      /NULCOM/ NCOUNT,USRNUL
C
      DATA     PI /3.1415926535/                         !Pi
C
C  branch according to function 
C
C  SQRT
      IF (CFUNC.EQ.'SQRT ') THEN
         IF (AVAL.LT.0.D0) THEN
            CVAL = DBLE(USRNUL)
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = DSQRT(AVAL)
         ENDIF
C
C  LN
      ELSE IF (CFUNC.EQ.'LN   ') THEN
         IF (AVAL.LE.0.D0) THEN
            CVAL = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = DLOG(AVAL)
         ENDIF
C
C  LOG10 or LOG
      ELSE IF (CFUNC(1:3).EQ.'LOG') THEN
         IF (AVAL.LE.0.D0) THEN
            CVAL = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = DLOG10(AVAL)
         ENDIF
C
C  EXP
      ELSE IF (CFUNC.EQ.'EXP  ') THEN
         CVAL = DEXP(AVAL)
C
C  EXP10
      ELSE IF (CFUNC.EQ.'EXP10') THEN
         CVAL = 10.D0 ** AVAL
C
C  SIN
      ELSE IF (CFUNC.EQ.'SIN  ') THEN
         TEMP = (AVAL / 180.D0) * PI                !degrees -> radians
         CVAL = DSIN(TEMP)
C
C  COS
      ELSE IF (CFUNC.EQ.'COS  ') THEN
         TEMP = (AVAL / 180.D0) * PI                !degrees -> radians
         CVAL = DCOS(TEMP)
C
C  TAN
      ELSE IF (CFUNC.EQ.'TAN  ') THEN
         TEMP = (AVAL / 180.D0) * PI                !degrees -> radians
         CVAL = DTAN(TEMP)
C
C  ASIN
      ELSE IF (CFUNC.EQ.'ASIN ') THEN
         TEMP = DASIN(AVAL)
         CVAL = (TEMP * 180.D0) / PI                !radians -> degrees
C
C  ACOS
      ELSE IF (CFUNC.EQ.'ACOS ') THEN
         TEMP = DACOS(AVAL)
         CVAL = (TEMP * 180.D0) / PI                !radians -> degrees
C
C  ATAN
      ELSE IF (CFUNC.EQ.'ATAN ') THEN
         TEMP = DATAN(AVAL)
         CVAL = (TEMP * 180.D0) / PI                !radians -> degrees
C
C  NINT
      ELSE IF (CFUNC.EQ.'INT  ') THEN
         CVAL = DBLE(NINT(AVAL))
C
C  ABS
      ELSE IF (CFUNC.EQ.'ABS  ') THEN
         CVAL = DABS(AVAL)
      ENDIF
C
      RETURN          
      END

      SUBROUTINE DF2CC(CFUNC,AVAL,CVAL)
C
      IMPLICIT NONE
C
      DOUBLE PRECISION   AVAL(2),CVAL,PI
      DOUBLE PRECISION   TEMP
C
      CHARACTER*5      CFUNC
C
      DATA     PI /3.1415926535/                         !Pi
C
C  branch according to function
C
C  function ATAN2
      IF (CFUNC.EQ.'ATAN2') THEN
         TEMP = DATAN2(AVAL(1),AVAL(2))
         CVAL = (TEMP * 180.D0) / PI                !radians -> degrees
C
C  get minimum
      ELSE IF (CFUNC.EQ.'MIN  ') THEN
         CVAL = MIN(AVAL(1),AVAL(2))
C
C  get maximum
      ELSE IF (CFUNC.EQ.'MAX  ') THEN
         CVAL = MAX(AVAL(1),AVAL(2))
C
C  get remainder
      ELSE IF (CFUNC.EQ.'MOD  ') THEN
         CVAL = DMOD(AVAL(1),AVAL(2))
      ENDIF
C
      RETURN
      END

      SUBROUTINE DOPFW(OPCODE,A,B,C,APIX,BPIX,CPIX,
     +                 NPIXA,NPIXB,NPIXC)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine OPFFW            version 2.00      880202
C             FN2FFW                   2.00      880202
C  K. Banse                        ESO - Garching
C  2.10		890915
C
C.KEYWORDS
C  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 OPFFW(OPCODE,A,B,C,APIX,BPIX,CPIX,NPIXA,NPIXB,NPIXC)
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  APIX:      I*4 array      start + end pixels for A
C  BPIX:      I*4 array      start + end pixels for B
C  CPIX:      I*4 array      start + end pixels for C
C  NPIXA:      I*4            no. of pixels per line in A
C  NPIXB:      I*4            no. of pixels per line in B
C  NPIXC:      I*4            no. of pixels per line in C
C
C  output par:
C  C:            R*8 array      result array
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      N,NN,NNN,NCOUNT
      INTEGER      NX,NY,NZ,OFFA,OFFB,OFFC
      INTEGER      XOFFA,XOFFB,XOFFC
      INTEGER      YOFFA,YOFFB,YOFFC
      INTEGER      ZOFFA,ZOFFB,ZOFFC
      INTEGER      APIX(3,2),BPIX(3,2),CPIX(3,2)
      INTEGER      NPIXA(*),NPIXB(*),NPIXC(*)
C
      CHARACTER*(*) OPCODE
C 
      REAL   USRNUL
C
      DOUBLE PRECISION     A(*),B(*),C(*)
      DOUBLE PRECISION     DEPSP,DEPSN
C
      COMMON      /NULCOM/ NCOUNT,USRNUL
C
      DATA DEPSP /10.D-33/
      DATA DEPSN /-10.D-33/
C
C  init
      NX = APIX(1,2) - APIX(1,1) + 1
      NY = APIX(2,2) - APIX(2,1) + 1
      NZ = APIX(3,2) - APIX(3,1) + 1
C
      XOFFA = APIX(1,1) - 1
      XOFFC = CPIX(1,1) - 1
      YOFFA = (APIX(2,1)-1)*NPIXA(1)
      YOFFC = (CPIX(2,1)-1)*NPIXC(1)
      ZOFFA = (APIX(3,1)-1)*NPIXA(1)*NPIXA(2)
      ZOFFC = (CPIX(3,1)-1)*NPIXC(1)*NPIXC(2)
C
C  copy 1. frame into result frame
      IF (OPCODE(3:3).EQ.'=') THEN
         DO 5500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 5300, NN = 1,NY
               DO 5200, N=1,NX
                  C(OFFC+N) = A(OFFA+N)
5200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFC = OFFC + NPIXC(1)
5300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
5500     CONTINUE
         RETURN
      ENDIF
C
C 
      XOFFB = BPIX(1,1) - 1            	!all other opcodes have 3 operands
      YOFFB = (BPIX(2,1)-1)*NPIXB(1)
      ZOFFB = (BPIX(3,1)-1)*NPIXB(1)*NPIXB(2)
C
C  branch according to opcode
C
C  add frames
      IF (OPCODE(3:3).EQ.'+') THEN
         DO 1500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 1300, NN = 1,NY
               DO 1100, N=1,NX
                  C(OFFC+N) = A(OFFA+N) + B(OFFB+N)
1100           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
1300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
1500     CONTINUE
C
C  subtract 2. frame from 1. frame
      ELSE IF (OPCODE(3:3).EQ.'-') THEN
         DO 2500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 2300, NN = 1,NY
               DO 2200, N=1,NX
                  C(OFFC+N) = A(OFFA+N) - B(OFFB+N)
2200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
2300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
2500     CONTINUE
C
C  multiply frames
      ELSE IF (OPCODE(3:3).EQ.'*') THEN
         DO 3500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 3300, NN = 1,NY
               DO 3200, N=1,NX
                  C(OFFC+N) = A(OFFA+N) * B(OFFB+N)
3200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
3300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
3500     CONTINUE
C
C  divide 1. frame by 2. frame
      ELSE IF (OPCODE(3:3).EQ.'/') THEN
C 
         DO 4500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 4300, NN = 1,NY
               DO 4200, N=1,NX
                  IF ((B(OFFB+N) .LT. DEPSP) .AND.
     +               (B(OFFB+N) .GT. DEPSN)) THEN
                     C(OFFC+N) = DBLE(USRNUL)
                     NCOUNT = NCOUNT + 1
                  ELSE
                     C(OFFC+N) = A(OFFA+N) / B(OFFB+N)
                  ENDIF
4200           CONTINUE
C 
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
4300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
4500     CONTINUE
C
C  compute 1. frame ** 2. frame
      ELSE IF (OPCODE(3:3).EQ.'`') THEN
         DO 6500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 6300, NN = 1,NY
               DO 6200, N=1,NX
                  C(OFFC+N) = A(OFFA+N) ** B(OFFB+N)
6200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
6300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
6500     CONTINUE
      ENDIF
C 
      RETURN
      END

      SUBROUTINE DF2FFW(CFUNC,A,B,C,APIX,BPIX,CPIX,
     +                  NPIXA,NPIXB,NPIXC)
C
      IMPLICIT NONE
C
      INTEGER      N,NN,NNN,NCOUNT
      INTEGER      NX,NY,NZ,OFFA,OFFB,OFFC
      INTEGER      XOFFA,XOFFB,XOFFC
      INTEGER      YOFFA,YOFFB,YOFFC
      INTEGER      ZOFFA,ZOFFB,ZOFFC
      INTEGER      APIX(3,2),BPIX(3,2),CPIX(3,2)
      INTEGER      NPIXA(*),NPIXB(*),NPIXC(*)
C
      CHARACTER*5   CFUNC
C
      DOUBLE PRECISION   A(*),B(*),C(*)
      DOUBLE PRECISION   FACT,RA,RB
C 
      REAL        USRNUL
C
      COMMON      /NULCOM/ NCOUNT,USRNUL
C
      DATA   FACT  /0.0174532925/
C
C  init
      NX = APIX(1,2) - APIX(1,1) + 1
      NY = APIX(2,2) - APIX(2,1) + 1
      NZ = APIX(3,2) - APIX(3,1) + 1
C
      XOFFA = APIX(1,1) - 1
      XOFFB = BPIX(1,1) - 1
      XOFFC = CPIX(1,1) - 1
      YOFFA = (APIX(2,1)-1)*NPIXA(1)
      YOFFB = (BPIX(2,1)-1)*NPIXB(1)
      YOFFC = (CPIX(2,1)-1)*NPIXC(1)
      ZOFFA = (APIX(3,1)-1)*NPIXA(1)*NPIXA(2)
      ZOFFB = (BPIX(3,1)-1)*NPIXB(1)*NPIXB(2)
      ZOFFC = (CPIX(3,1)-1)*NPIXC(1)*NPIXC(2)
C
C  branch according to function
C
C  function ATAN2
      IF (CFUNC.EQ.'ATAN2') THEN
         DO 1500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 1300, NN = 1,NY
               DO 1200, N=1,NX
                  RA = A(OFFA+N) * FACT
                  RB = B(OFFB+N) * FACT
                  C(OFFC+N) = DATAN2(RA,RB)
1200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
1300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
1500     CONTINUE
C
C  minimum of 2 frames
      ELSE IF (CFUNC.EQ.'MIN  ') THEN
         DO 2500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 2300, NN = 1,NY
               DO 2200, N=1,NX
                  C(OFFC+N) = MIN(A(OFFA+N),B(OFFB+N))
2200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
2300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
2500     CONTINUE
C
C  maximum of 2 frames
      ELSE IF (CFUNC.EQ.'MAX  ') THEN
         DO 3500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 3300, NN = 1,NY
               DO 3200, N=1,NX
                  C(OFFC+N) = MAX(A(OFFA+N),B(OFFB+N))
3200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
3300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
3500     CONTINUE
C
C  remainder of 1.frame/2. frames
      ELSE IF (CFUNC.EQ.'MOD  ') THEN
         DO 4500, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 4300, NN = 1,NY
               DO 4200, N=1,NX
                  C(OFFC+N) = DMOD(A(OFFA+N),B(OFFB+N))
4200           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
4300        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
4500     CONTINUE
      ENDIF
C 
      RETURN
      END
