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 OPCC(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.FOR
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*4 array      1. input array
C  B:            R*4 array      2. input array
C  SCALR:      R*4            scalar value
C  CONSTA:      R*4            1. input scalar
C  CONSTB:      R*4            2. input scalar
C  NDIM:      I*4            size of arrays involved
C
C  output par:
C  C:            R*4 array      result array
C  CONSTC:      R*4            result scalar
C
C.VERSION
C 030710	last modif
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NCOUNT
C
      REAL     CONSTA,CONSTB,CONSTC
      REAL     USRNUL,EPSP,EPSN
C
      CHARACTER*(*)      OPCODE
C
      COMMON      /NULCOM/NCOUNT,USRNUL
C 
      DATA  EPSP /10.0E-33/, EPSN /-10.0E-33/
C
C  branch according to opcode
C
C  subtract 2. scalar from 1. scalar value
      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.EPSP) .AND.
     +       (CONSTB.GT.EPSN)) 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
C	on OSF/DEC Unix the usual FORTRAN statement:
C	CONSTC = CONSTA ** CONSTB 
C	does not work if CONSTA < 0. - so we have to do it in C ...
         CALL CPOWER(CONSTA,CONSTB,CONSTC)
      ELSE
         WRITE(*,10000)
      ENDIF
      RETURN
C
10000 FORMAT(' OPCC: We should not come here...')
C
      END

      SUBROUTINE FUN1C(CFUNC,AVAL,CVAL)
C
      IMPLICIT NONE
C
      INTEGER    NCOUNT
C
      REAL       AVAL,CVAL
      REAL       TEMP,USRNUL,PI
C
      CHARACTER*5  CFUNC
C
      COMMON      /NULCOM/NCOUNT,USRNUL
C
      DATA     PI /3.141593/                         !Pi
C
C  branch according to function 
C
C  SQRT
      IF (CFUNC.EQ.'SQRT ') THEN
         IF (AVAL.LT.0.) THEN
            CVAL = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = SQRT(AVAL)
         ENDIF
C
C  LN
      ELSE IF (CFUNC.EQ.'LN   ') THEN
         IF (AVAL.LE.0.) THEN
            CVAL = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = LOG(AVAL)
         ENDIF
C
C  LOG10 or LOG
      ELSE IF (CFUNC(1:3).EQ.'LOG') THEN
         IF (AVAL.LE.0.) THEN
            CVAL = USRNUL
            NCOUNT = NCOUNT + 1
         ELSE
            CVAL = LOG10(AVAL)
         ENDIF
C
C  EXP
      ELSE IF (CFUNC.EQ.'EXP  ') THEN
         CVAL = EXP(AVAL)
C
C  EXP10
      ELSE IF (CFUNC.EQ.'EXP10') THEN
         CVAL = 10. ** AVAL
C
C  SIN
      ELSE IF (CFUNC.EQ.'SIN  ') THEN
         TEMP = (AVAL / 180.) * PI                !degrees -> radians
         CVAL = SIN(TEMP)
C
C  COS
      ELSE IF (CFUNC.EQ.'COS  ') THEN
         TEMP = (AVAL / 180.) * PI                !degrees -> radians
         CVAL = COS(TEMP)
C
C  TAN
      ELSE IF (CFUNC.EQ.'TAN  ') THEN
         TEMP = (AVAL / 180.) * PI                !degrees -> radians
         CVAL = TAN(TEMP)
C
C  ASIN
      ELSE IF (CFUNC.EQ.'ASIN ') THEN
         TEMP = ASIN(AVAL)
         CVAL = (TEMP * 180.) / PI                !radians -> degrees
C
C  ACOS
      ELSE IF (CFUNC.EQ.'ACOS ') THEN
         TEMP = ACOS(AVAL)
         CVAL = (TEMP * 180.) / PI                !radians -> degrees
C
C  ATAN
      ELSE IF (CFUNC.EQ.'ATAN ') THEN
         TEMP = ATAN(AVAL)
         CVAL = (TEMP * 180.) / PI                !radians -> degrees
C
C  NINT
      ELSE IF (CFUNC.EQ.'INT  ') THEN
         CVAL = FLOAT(NINT(AVAL))
C
C  ABS
      ELSE IF (CFUNC.EQ.'ABS  ') THEN
         CVAL = ABS(AVAL)
C
      ELSE
         WRITE(*,10000)
      ENDIF
      RETURN          
C 
10000 FORMAT(' FUN1C: We should not come here...')
      END

      SUBROUTINE FUN2CC(CFUNC,AVAL,CVAL)
C
      IMPLICIT NONE
C
      REAL          AVAL(2),CVAL,PI
      REAL          TEMP
C
      CHARACTER*5      CFUNC
C
      DATA     PI /3.141593/                         !Pi
C
C  branch according to function
C
C  function ATAN2
      IF (CFUNC.EQ.'ATAN2') THEN
         TEMP = ATAN2(AVAL(1),AVAL(2))
         CVAL = (TEMP * 180.) / 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 = MOD(AVAL(1),AVAL(2))
C
      ELSE
         WRITE(*,10000)
      ENDIF
      RETURN
C
10000 FORMAT(' FUN2CC: We should not come here...')
C
      END

      SUBROUTINE OPFFW(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
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*4 array      1. input array
C  B:            R*4 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*4 array      result array
C
C.VERSION
C 030710	last modif
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     A(*),B(*),C(*)
      REAL     USRNUL,EPSP,EPSN,RR
C
      COMMON      /NULCOM/NCOUNT,USRNUL
C
      DATA  EPSP /10.0E-33/, EPSN /-10.0E-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 4800, NNN=1,NZ
            OFFA = ZOFFA + YOFFA + XOFFA
            OFFB = ZOFFB + YOFFB + XOFFB
            OFFC = ZOFFC + YOFFC + XOFFC
C
            DO 4700, NN = 1,NY
               DO 4600, N=1,NX
                  RR = B(OFFB+N)
                  IF ((RR .LT. EPSP) .AND. (RR .GT. EPSN)) THEN
                     C(OFFC+N) = USRNUL
                     NCOUNT = NCOUNT + 1
                  ELSE
                     C(OFFC+N) = A(OFFA+N) / RR
                  ENDIF
4600           CONTINUE
               OFFA = OFFA + NPIXA(1)
               OFFB = OFFB + NPIXB(1)
               OFFC = OFFC + NPIXC(1)
4700        CONTINUE
C
            ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2)
            ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2)
            ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2)
4800     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
C 
      ELSE
         WRITE(*,10000)
      ENDIF
      RETURN
C
10000 FORMAT(' OPFFW: We should not come here...')
      END

      SUBROUTINE FN2FFW(CFUNC,A,B,C,APIX,BPIX,CPIX,
     +                  NPIXA,NPIXB,NPIXC)
C
      IMPLICIT NONE
C
      INTEGER      N,NN,NNN
      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
      REAL          A(*),B(*),C(*),FACT,RA,RB
C
      DATA   FACT  /0.0174533/
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) = ATAN2(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) = MOD(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
C 
      ELSE
         WRITE(*,10000)
      ENDIF
      RETURN
C
10000 FORMAT(' FN2FFW: We should not come here...')
      END
C 

      SUBROUTINE DSPNUL(NCC)
C 
      INTEGER     NCC
      INTEGER     NCOUNT,STAT
C 
      CHARACTER   OUTPUT*80
C 
      REAL        USRNUL
      COMMON      /NULCOM/NCOUNT,USRNUL
C      
      IF (NCC.GT.1) THEN
         WRITE(OUTPUT,10000) NCC,USRNUL
      ELSE
         WRITE(OUTPUT,10001) USRNUL
      ENDIF
      CALL STTPUT(OUTPUT,STAT)
      RETURN
C      
C  format
C 
10000 FORMAT(I7,' undefined pixels ... set to "null value" = ',G15.7)
10001 FORMAT('1 undefined pixel ... set to "null value" = ',G15.7)
C      
      END
