C @(#)brthmz.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:31
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 
     +COMPUB(RESFMT,EXPRSS,ATOM,APNTRS,UPDA,FRAMEC,CUTS)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine COMPUB         version 4.50      880920
C  K. Banse                  ESO - Garching
C  				     4.60      890426
C
C.KEYWORDS
C  arithmetic operations, bulk data frames
C
C.PURPOSE
C  evaluate an arithmetic expression involving frames, constants and functions
C  and store the result into a data frame or key OUTPUTR(1), if used as pocket calculator
C
C.ALGORITHM
C  "clean" the expression by replacing all frame names by F, all constants by C
C  and all functions by P, convert it to polish (postfix) notation and evaluate it piecewise
C  intermediate results are stored in OUTPUTR(1) or frame midtempn.bdf
C  establish special condition handler for arithmetic traps...
C  copy all descriptors of first frame operand to result frame.
C
C.INPUT/OUTPUT
C      
C  call as  COMPUB(RESFMT,EXPRSS,ATOM,APNTRS,UPDA,FRAMEC,CUTS)
C 
C  input par:
C  RESFMT:	I*4           data format for result frame 
C  EXPRSS:      char.exp.     arithmetic expression in polish postfix notation
C  ATOM:        char. array   holds the operands
C  APNTRS:      I*4 array     holds the pointers to ATOM
C  FRAMEC:      char. exp.    name of result frame
C  UPDA:	I*4           if = 1, result frame is just modified,
C                             else a new result frame is created
C      
C  output par:
C  CUTS:        R*4 array     min, max of result frame
C      
C-------------------------------------------------------------------------
C
      IMPLICIT NONE
C
      CHARACTER*(*)    EXPRSS,ATOM(*),FRAMEC
      CHARACTER*80     FRAMEA,FRAMEB,FRAME
      CHARACTER*1      DUM(24)
      CHARACTER        CUNITA*64,CUNITB*64,IDENTA*72,IDENTB*72
      CHARACTER*60     OPERA,OPERB,OPERC
      CHARACTER        WORK(2)*50
      CHARACTER        OPERAT*4,LINE*80,DUMMY*20
      CHARACTER        ERROR1*40,ERROR2*30,ERROR4*30
C      
      DOUBLE PRECISION STEPA(3),STEPB(3),STEPC(3)
      DOUBLE PRECISION STARTA(3),STARTB(3),STARTC(3)
      DOUBLE PRECISION BEGIN(3),END(3),DIF,DDUM(1)
C 
      REAL         CNSTAB(2),CONST,EPS1,CUTS(2),RDUM(1)
C      
      INTEGER      RESFMT,APNTRS(*),UPDA,IMNOC
      INTEGER      APIX(3,2),BPIX(3,2),CPIX(3,2)
      INTEGER      NPIXA(3),NPIXB(3),NPIXC(3)
      INTEGER      IDUM(1),IDUMMY,LL
      INTEGER      MAXDIM,N,NAXISA,NAXISB,NAXISC
      INTEGER      NBRA,P1,P2,P3
      INTEGER*8    PNTRA,PNTRB,PNTRC
      INTEGER      PP,SIZE,STAT
      INTEGER      IMNOA,IMNOB
      INTEGER      UNI(1),NULO,MADRID(1)
      INTEGER      APIXDF,BPIXDF,CPIXDF,KK
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C 
      COMMON      /VMR/ MADRID
C      
      DATA      
     +ERROR1 /'Operands do not match in stepsize... '/
      DATA      ERROR2      /'Operands do not overlap...'/
      DATA      ERROR4      /'Too many operands...'/
c 
      DATA      APIX /6*1/, BPIX /6*1/, CPIX /6*1/
      DATA      DUM  /'a','b','c','d','e','f','g','h','i','j',
     +                'k','l','m','n','o','p','q','r','s','t',
     +                'u','v','w','x'/
      DATA      IDENTA /' '/,  CUNITA /' '/
      DATA      IDENTB /' '/,  CUNITB /' '/
      DATA      NPIXA  /3*1/,    NPIXB /3*1/,     NPIXC /3*1/
      DATA      STARTA /3*0.D0/, STARTB /3*0.D0/, STARTC /3*0.D0/
      DATA      STEPA  /3*1.D0/, STEPB  /3*1.D0/, STEPC  /3*1.D0/
      DATA      DUMMY  /'midtemp   .bdf '/
C 
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C      
C  initialize
      IDUMMY = 0
      WORK(1) = EXPRSS
      CALL STKRDC('MID$SESS',1,11,2,LL,DUMMY(8:),UNI,NULO,STAT)
C      
C  extract basic operations
1000  CALL EXPRDC(WORK(1),WORK(2),OPERAT,PP)
C      
C  extract operands
      P1 = APNTRS(PP)
      P2 = APNTRS(PP+1)
      OPERA = ATOM(P1)
      OPERB = ATOM(P2)
      IF (OPERAT(1:1).EQ.'Q') THEN                  !treat 2-arg functions...
         P3 = APNTRS(PP+2)
         OPERC = ATOM(P3)
      ENDIF
C      
C  find out what kind of OPERA to do
      IF (INDEX(OPERAT,'C').GT.0) THEN
         N = INDEX(OPERAT,'F')                  !also F involved...?
         IF (N.GT.0) GOTO 1600                  !yes. 
      ELSE
         IF (OPERAT(1:1).EQ.'P') THEN
            GOTO 1600                              !handle PF)
         ELSE
            GOTO 2000                              !handle FF or QFF)
         ENDIF
      ENDIF
C      
C  ***
C  only constants involved, do it right now
C  ***
C      
      IF (OPERAT(1:2).EQ.'CC') THEN
         LL = INDEX(OPERA,' ') - 1
         IF (LL.LE.0) LL = 60
         LINE(1:) = OPERA(1:LL)//','//OPERB
         CALL GENCNV(LINE(1:61),2,2,IDUM,CNSTAB,DDUM,LL)
         IF (LL.LE.1) GOTO 9990                              !we need 2 values...
         CALL OPCC(OPERAT,CNSTAB(1),CNSTAB(2),CONST)
      ELSE
         IF (OPERAT(1:1).EQ.'P') THEN                  !1-arg functions
            CALL GENCNV(OPERB,2,1,IDUM,CNSTAB,DDUM,LL)
            IF (LL.LE.0) GOTO 9990
            CALL FUN1C(OPERA(1:5),CNSTAB,CONST)
         ELSE                                          !2-arg functions
            LINE(1:41) = OPERB//','//OPERC
            CALL GENCNV(LINE(1:61),2,2,IDUM,CNSTAB,DDUM,LL)
            IF (LL.LE.0) GOTO 9990
            CALL FUN2CC(OPERA(1:5),CNSTAB,CONST)
         ENDIF
      ENDIF
C
C  put resulting constant back into relevant ATOM + goto loopend
      IF (WORK(2)(2:2).NE.' ') THEN
         WRITE(ATOM(P1),10000) CONST
         GOTO 5000
      ELSE
         GOTO 9000
      ENDIF
C
C  ***
C  one operand is a file
C  ***
C      
1600  IF (OPERAT(1:2).EQ.'FC') THEN
         FRAME = OPERA
         LINE(1:20) = OPERB
         NBRA = 1                                    !no function
         GOTO 1700
      ENDIF
C 
      IF (OPERAT(1:2).EQ.'CF') THEN
         FRAME = OPERB
         LINE(1:20) = OPERA
         NBRA = 1                                    !no function
         GOTO 1700
      ENDIF
C 
      IF (OPERAT(1:2).EQ.'PF') THEN
         FRAME = OPERB
         NBRA = 2                                    !1-arg function
         GOTO 1800
      ENDIF
      NBRA = 3                                    !2-arg function
C 
      IF (OPERAT(1:2).NE.'QC') THEN
         FRAME = OPERB
         LINE(1:20) = OPERC
      ELSE
         FRAME = OPERC
         LINE(1:20) = OPERB
      ENDIF
C      
1700  CALL GENCNV(LINE(1:20),2,1,IDUM,RDUM,DDUM,LL)
      CONST = RDUM(1)
      IF (LL.LE.0) GOTO 9990                     !problems with conversion...
C      
1800  CALL CLNFRA(FRAME,FRAMEA,0)                !take care of & and # ...
C      
C  now map the result frame
      IF ( (WORK(2)(2:2).EQ.' ') .AND.
     +     (FRAMEA.EQ.FRAMEC) ) THEN
         CALL STIGET(FRAMEA,RESFMT,F_IO_MODE,F_IMA_TYPE,
     +               3,NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         SIZE = 1
C 
         DO 1805, N=1,NAXISC
            SIZE = SIZE * NPIXC(N)
1805     CONTINUE
         PNTRC = PNTRA
         IMNOC = IMNOA
      ELSE
C      
C  result frame different from input frames
         CALL STIGET(FRAMEA,RESFMT,F_I_MODE,F_IMA_TYPE,
     +               3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         NAXISC = NAXISA
         SIZE = 1
C 
         DO 1810, N=1,NAXISA
            SIZE = SIZE * NPIXA(N)
            NPIXC(N) = NPIXA(N)
            STARTC(N) = STARTA(N)
            STEPC(N) = STEPA(N)
1810     CONTINUE
C      
C  last operation ?
         IF (WORK(2)(2:2).EQ.' ') THEN
            IF (UPDA.EQ.1) THEN         
               CALL STIGET(FRAMEC,RESFMT,F_IO_MODE,F_IMA_TYPE,
     +                     3,NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +                     CUNITA,PNTRC,IMNOC,STAT)
C                                            we have to recalculate SIZE
               SIZE = 1                       
               DO 1815, N=1,NAXISC
                  SIZE = SIZE * NPIXC(N)
1815           CONTINUE
            ELSE
               CALL STIPUT(FRAMEC,RESFMT,F_O_MODE,F_IMA_TYPE,
     +                     NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +                     CUNITA,PNTRC,IMNOC,STAT)
            ENDIF
         ELSE                                      !use dummy result frame
C       
C  no. loop more
            IDUMMY = IDUMMY + 1                    !increment dummy file counter
            IF (IDUMMY.GT.24) THEN                 !check operand count
               CALL STETER(3,ERROR4)
            ELSE
               DUMMY(10:10) = DUM(IDUMMY)
            ENDIF
            CALL STIPUT(DUMMY,RESFMT,F_O_MODE,F_IMA_TYPE,
     +                  NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +                  CUNITA,PNTRC,IMNOC,STAT)
         ENDIF
      ENDIF
C
C  now do the actual operation
      GOTO (1850,1860,1870),NBRA                       
C 
1850  CALL OPFC(OPERAT,MADRID(PNTRA),CONST,MADRID(PNTRC),SIZE)
      GOTO 1900
C 
1860  CALL FUN1F(OPERA(1:5),MADRID(PNTRA),MADRID(PNTRC),SIZE)
      GOTO 1900
C 
1870  CALL FUN2FC(OPERA(1:5),MADRID(PNTRA),CONST,MADRID(PNTRC),SIZE)
C      
1900  IF (WORK(2)(2:2).NE.' ') THEN
         CALL STFCLO(IMNOA,STAT) 
         GOTO 4400                           
      ELSE
         GOTO 9000                           
      ENDIF
C
C  ***
C  both operands are files 
C  ***
C      
2000  IF (OPERAT(1:1).NE.'Q') THEN
         CALL CLNFRA(OPERA,FRAMEA,0)
         CALL CLNFRA(OPERB,FRAMEB,0)
      ELSE
         CALL CLNFRA(OPERB,FRAMEA,0)
         CALL CLNFRA(OPERC,FRAMEB,0)
      ENDIF
C      
C  get first input frame
      IF ( (WORK(2)(2:2).EQ.' ') .AND.
     +     (FRAMEA.EQ.FRAMEC) ) THEN
         CALL STIGET(FRAMEA,RESFMT,F_IO_MODE,F_IMA_TYPE,
     +               3,NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)      !map input/result frame
         NAXISA = NAXISC
         SIZE = 1
C 
         DO 2050, N=1,NAXISC
            SIZE = SIZE * NPIXC(N)
            NPIXA(N) = NPIXC(N)
            STARTA(N) = STARTC(N)
            STEPA(N) = STEPC(N)
2050     CONTINUE
         PNTRC = PNTRA
         IMNOC = IMNOA
      ELSE
         CALL STIGET(FRAMEA,RESFMT,F_I_MODE,F_IMA_TYPE,
     +               3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
      ENDIF
C      
C  handle 2. input frame
      IF (FRAMEA.EQ.FRAMEB) THEN                        !check, if both operands are the same
         PNTRB = PNTRA
         IMNOB = IMNOA
         DO 2100, N=1,NAXISA
            NPIXB(N) = NPIXA(N)
            STARTB(N) = STARTA(N)
            STEPB(N) = STEPA(N)
2100     CONTINUE
         NAXISB = NAXISA
      ELSE
         CALL STIGET(FRAMEB,RESFMT,F_I_MODE,F_IMA_TYPE,
     +               3,NAXISB,NPIXB,STARTB,STEPB,IDENTB,
     +               CUNITB,PNTRB,IMNOB,STAT)
      ENDIF
C
C  see, if stepsizes and origins fit + frames overlap 
      IF (NAXISB.GT.NAXISA) THEN
         MAXDIM = NAXISB                        !take maximum
      ELSE
         MAXDIM = NAXISA
      ENDIF
C 
      DO 2150, N=1,MAXDIM
         EPS1 = 0.001 * ABS(STEPA(N))                  !take 0.1% of stepsize as epsilon
         IF (ABS(STEPA(N)-STEPB(N)).GT.EPS1) 
     +      CALL STETER(1,ERROR1)
         CALL OVRLAP(STARTA(N),STEPA(N),NPIXA(N),STARTB(N),
     +               STEPB(N),NPIXB(N),BEGIN(N),END(N),STAT)
         IF (STAT.NE.0) CALL STETER(3,ERROR2)        !if STAT = 1, no overlap...
2150  CONTINUE
C 
C  create new result frame with dimension = intersection of input frames
      IF ( (WORK(2)(2:2).EQ.' ') .AND.
     +     (FRAMEA.EQ.FRAMEC) ) GOTO 3000            !we already have the result frame
C      
      IF (NAXISA.GT.NAXISB) THEN
         NAXISC = NAXISB
      ELSE
         NAXISC = NAXISA
      ENDIF
C 
      SIZE = 1
      DO 2200, N=1,NAXISC
         STARTC(N) = BEGIN(N)
         STEPC(N) = STEPA(N)
         NPIXC(N) = NINT( (END(N)-BEGIN(N)) / STEPC(N) ) + 1
         SIZE = SIZE * NPIXC(N)
2200  CONTINUE
C      
C  now map the resulting frame
      IF (WORK(2)(2:2).EQ.' ') THEN                !use real result frame
         IF (UPDA.EQ.1) THEN                      !test, if update or new file
            CALL STIGET(FRAMEC,RESFMT,F_IO_MODE,F_IMA_TYPE,
     +                  3,NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +                  CUNITA,PNTRC,IMNOC,STAT)
            SIZE = 1                              !we have to recalculate SIZE
            DO 2250, N=1,NAXISC
               SIZE = SIZE * NPIXC(N)
2250        CONTINUE
         ELSE
            CALL STIPUT(FRAMEC,RESFMT,F_O_MODE,F_IMA_TYPE,
     +                  NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +                  CUNITA,PNTRC,IMNOC,STAT)
         ENDIF
C      
      ELSE                                         !use dummy result frame
         IDUMMY = IDUMMY + 1                       !increment dummy file counter
         IF (IDUMMY.GT.24) THEN                    !check operand count
            CALL STETER(2,ERROR4)
         ELSE
            DUMMY(10:10) = DUM(IDUMMY)
         ENDIF
         CALL STIPUT(DUMMY,RESFMT,F_O_MODE,F_IMA_TYPE,
     +               NAXISC,NPIXC,STARTC,STEPC,IDENTA,
     +               CUNITA,PNTRC,IMNOC,STAT)
      ENDIF
C 
C  convert start + end of overlap region into pixel no.'s
3000  DO 3100, N=1,MAXDIM
         DIF = (BEGIN(N)-STARTA(N))/STEPA(N)
         APIX(N,1) = NINT(DIF) + 1
         DIF = (END(N)-STARTA(N))/STEPA(N)
         APIX(N,2) = NINT(DIF) + 1
         APIXDF = APIX(N,2)-APIX(N,1)
         DIF = (BEGIN(N)-STARTB(N))/STEPB(N)
         BPIX(N,1) = NINT(DIF) + 1
         DIF = (END(N)-STARTB(N))/STEPB(N)
         BPIX(N,2) = NINT(DIF) + 1
         BPIXDF = BPIX(N,2)-BPIX(N,1)
         DIF = (BEGIN(N)-STARTC(N))/STEPC(N)
         CPIX(N,1) = NINT(DIF) + 1
         DIF = (END(N)-STARTC(N))/STEPC(N)
         CPIX(N,2) = NINT(DIF) + 1
         CPIXDF = CPIX(N,2)-CPIX(N,1)
         KK = MIN(APIXDF,BPIXDF,CPIXDF)            !take smallest size
         APIX(N,2) = APIX(N,1) + KK
         BPIX(N,2) = BPIX(N,1) + KK
         CPIX(N,2) = CPIX(N,1) + KK
3100  CONTINUE
C 
C  now do the actual operation
      IF (OPERAT(1:1).NE.'Q') THEN
         CALL OPFFW(OPERAT,MADRID(PNTRA),MADRID(PNTRB),
     +          MADRID(PNTRC),APIX,BPIX,CPIX,NPIXA,NPIXB,NPIXC)
      ELSE
         CALL FN2FFW(OPERA(1:5),MADRID(PNTRA),MADRID(PNTRB),
     +          MADRID(PNTRC),APIX,BPIX,CPIX,NPIXA,NPIXB,NPIXC)
      ENDIF
      IF (WORK(2)(2:2).EQ.' ') GOTO 9000
C      
      CALL STFCLO(IMNOA,STAT)      !release input frames
      CALL STFCLO(IMNOB,STAT) 
C
C  put resulting frame back into relevant ATOM, if we are not finished yet
4400  ATOM(P1) = DUMMY
      CALL STFCLO(IMNOC,STAT)      !release dummy result frame
C      
C  loopend for all basic operations
5000  IF (OPERAT(1:1).NE.'Q') THEN
         NBRA = 2
      ELSE
         NBRA = 3
      ENDIF
C 
      DO 5050, N=PP+1,45                        !update pointers for ATOM
         APNTRS(N) = APNTRS(N+NBRA)
5050  CONTINUE
C 
      WORK(1) = WORK(2)
      GOTO 1000                        !get next operation
C  
C  ***
C  we're done. test, if result goes to frame or constant
C  ***
C      
C  calculate new dynamic range of result frame
9000  CALL MYMX(MADRID(PNTRC),SIZE,CUTS)
      CALL STFCLO(IMNOC,STAT)                 !make sure result frame is closed
C      
C  delete any dummy frames
      IF (IDUMMY.GT.0) THEN
         DO 9100, N=1,IDUMMY
            DUMMY(10:10) = DUM(N)
            CALL STFDEL(DUMMY,STAT)
9100     CONTINUE
      ENDIF
C      
C  That's it folks
      RETURN
C      
C  error with conversion of ASCII to number
9990  CALL STETER(10,'conversion error of ASCII -> number ...')
C
C  formats...
10000 FORMAT(G15.7)
      END