C @(#)arthmd.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:30
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 ARTHMD(ZDEFA,ZRESFR,ZLINE)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine ARTHMD			951011
C  K. Banse            		      	ESO - Garching
C
C.KEYWORDS
C  double precision, 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), 
C  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  establish special condition handler for arithmetic traps...
C  copy all descriptors of first frame operand to result frame.
C
C.INPUT/OUTPUT
C  call as ARTHMD(ZDEFA,ZRESFR,ZLINE)
C  
C  input par:
C  ZDEFA:	char.		2 char. default string
C				(1) Y or N for result frame or not
C				(2) P, I, X for COMPUTE/PIXEL, /IMAGE, /PLANE
C  ZRESFR:	char.		name of result frame
C  ZLINE:	char.		char.string with arithmetic expression
C	
C.VERSIONS
C	
C 000621
C-------------------------------------------------------------------------
C
      IMPLICIT NONE
C
C  important parameter: WORKLM = 45
C
      CHARACTER*(*)  ZRESFR,ZLINE
      CHARACTER      ZDEFA*2
      CHARACTER      INFRAM(23)*80                   !MAXFRAM = WORKLM/2 + 1
      CHARACTER      ATOM(23)*80,OPER(23)*5          !size = MAXFRAM
      CHARACTER      RESFRA*80,DSCFRA*80
      CHARACTER      CUNIT*80,IDENT*72
      CHARACTER      WORK*50,WORK1*50,OUTPU*80,WORKY*1
      CHARACTER      ERROR3*50,ERROR4*30,ERROR5*30
C      
      REAL           RCUTS(4),RDUM(1) 
C 
      DOUBLE PRECISION STEP(6),START(6),STEPA(6),STARTA(6)
      DOUBLE PRECISION DATTA(23),EPS,DIF,DDUM(1)
      DOUBLE PRECISION CONST,CUTS(4)
C      
      INTEGER      LATOM(23),FELEM(23),WORKLM        !size = MAXFRAM
      INTEGER      APNTRS(48)                        !size = WORKLM + 3
      INTEGER*8    PNTR(27)                          !size = MAXFRAM + 4
      INTEGER*8    RPNTR
      INTEGER      PIDX(23),PPP(23),IMNO(23)         !size = MAXFRAM
      INTEGER      FRACNT,NPIX(6),NAXISA,NPIXA(6)
      INTEGER      COUNT,IAV,K,M,LL
      INTEGER      N,NAXIS,TOTAL
      INTEGER      SIZE,INFO(5),DATFMT
      INTEGER      RFELEM,OPCNT
      INTEGER      MAPSIZ,STAT
      INTEGER      TRSIZE,NNULL,UNIT(1),INPU(1)
      INTEGER      RIMNO,WIMNO,MADRID(1)
      INTEGER      UPDA,DSCORG
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C      
      COMMON      /IPOOL/  FRACNT
      COMMON      /VMR/ MADRID
C 
      SAVE
C 
      DATA      WORKLM      /45/
      DATA      CUNIT       /' '/,  IDENT       /' '/
      DATA      ERROR3
     +          /'wrong syntax in arithmetic expression ...'/
      DATA      ERROR4      /'too many operands...'/
      DATA      ERROR5      /'invalid plane no...'/
C      
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C      
C  init variables
      CUTS(1) = 0.
      CUTS(2) = 0.
      CUTS(3) = + 999999.9
      CUTS(4) = - 999999.9
      RCUTS(1) = 0.
      RCUTS(2) = 0.
      RCUTS(3) = CUTS(3)
      RCUTS(4) = CUTS(4)
      RIMNO = -1
      RFELEM = 1
C 
      DO 120, N=1,6
         NPIX(N) = 1
         START(N) = 0.D0
         STEP(N) = 1.D0
120   CONTINUE
C 
      DO 130, N=1,23
         DATTA(N) = 0.D0
         FELEM(N) = 1
130   CONTINUE
      DO 140, N=1,48
         APNTRS(N) = -1
140   CONTINUE
C 
C  reduce expression... (clean up + fill ATOM, LATOM)
200   SIZE = 80                                !length of ATOM(i) strings
      CALL EXPCLE(ZLINE,WORK,SIZE,COUNT,ATOM,LATOM)      
      IF (COUNT.EQ.-1) THEN
         CALL STETER(1,ERROR3)     
      ELSE IF (COUNT.EQ.-2) THEN
         CALL STETER(5,'arithmetic expression too long...')
C      
C  take care of operand without operation   (extend F or C to F+C or C+C)
      ELSE IF (COUNT.EQ.1) THEN    
         ATOM(2) = '0.'    
         LATOM(2) = 2
         WORK(2:3) = '+C'
         COUNT = 2
      ENDIF
C      
C  convert cleaned algebraic expression to polish postfix
      ZLINE(1:) = ' '
      CALL EXPPOL(WORK,ZLINE,STAT)                 
      IF (STAT.NE.0) CALL STETER(2,ERROR3)     
      LL = INDEX(ZLINE,' ') - 1
      IF (LL.GT.WORKLM) 
     +   CALL STETER(8,'Too many operands - reduce expression...')
C      
C  fill pointers for atoms
      DATFMT = D_R8_FORMAT                   !that's what we use currently
      K = 1
      DO 220, N=1,LL
         WORKY = ZLINE(N:N)
         IF ( (WORKY.EQ.'C').OR. (WORKY.EQ.'F').OR.
     +        (WORKY.EQ.'Q').OR. (WORKY.EQ.'P') ) THEN
            APNTRS(N) = K
            K = K + 1
         ENDIF
220   CONTINUE
C 
C  check for descriptor-origin-frame (other than 1. frame)
C  but be careful of frame[.. , ..:.. , ..] syntax
C 
      DSCFRA(1:) = ' '
      LL = 0
      K = INDEX(ZRESFR,']')
      IF (K.GT.1) THEN                          !abc[...],descframe possible
         K = K + 1
         IF (ZRESFR(K:K).EQ.',') LL = K
      ELSE
         LL = INDEX(ZRESFR,',')                 !abc,descframe possible
      ENDIF
      IF (LL.GT.1) THEN
         OUTPU(1:) = ZRESFR(LL+1:)
         CALL CLNFRA(OUTPU,DSCFRA,0)
         ZRESFR(LL:) = ' '
      ENDIF
C 
C  if COMPUTE/PLANE we have to work more
      IF (ZDEFA(2:2).EQ.'X') THEN
         M = INDEX(ZRESFR,'.P')                 !we know it's there...
         CALL GENCNV(ZRESFR(M+2:),1,1,INPU,RDUM,DDUM,STAT)
         IF (STAT.LT.1) CALL STETER(24,ERROR5)
         RFELEM = INPU(1)
         IF (RFELEM.LT.1) CALL STETER(24,ERROR5)
         ZRESFR(M:) = ' '
      ENDIF
      CALL CLNFRA(ZRESFR,RESFRA,0)
C      
C  find frames + constants
      LL = 0                        !counter in ATOM
      K = 0
C 
      DO 250, N=1,WORKLM
         WORKY = WORK(N:N)
         IF (WORKY.EQ.' ') GOTO 300
C 
         IF (WORKY.EQ.'F') THEN
            K = K + 1                   !not more than 23 frames possible
            LL = LL + 1
C 
            IF (ZDEFA(2:2).EQ.'X') THEN
               M = INDEX(ATOM(LL),'.P')
               IF (M.LT.2) M = INDEX(ATOM(LL),'.p')
               IF (M.LT.2) THEN                       !check, if it's p#
                  IF ((ATOM(LL)(1:1).EQ.'P').OR.
     +                (ATOM(LL)(1:1).EQ.'p')) THEN
                     CALL GENCNV(ATOM(LL)(2:),1,1,INPU,RDUM,DDUM,STAT)
                     IF (STAT.EQ.1) THEN
                        FELEM(K) = INPU(1)
                        INFRAM(K) = RESFRA
                        GOTO 240
                     ENDIF
                  ENDIF
               ELSE                                   !it's frame.p#
                  CALL GENCNV(ATOM(LL)(M+2:),1,1,INPU,RDUM,DDUM,STAT)
                  IF (STAT.LT.1) CALL STETER(24,ERROR5)
                  FELEM(K) = INPU(1)
                  ATOM(LL)(M:) = ' '
               ENDIF
            ENDIF
C 
            CALL CLNFRA(ATOM(LL),INFRAM(K),0)
240         PIDX(LL) = K                !we'll use PIDX instead of ATOM later
C      
         ELSE IF (WORKY.EQ.'C') THEN
            LL = LL + 1                 !we'll use DATTA instead of ATOM later
            CALL GENCNV(ATOM(LL),4,1,INPU,RDUM,DATTA(LL),STAT)
            IF (STAT.LE.0) CALL STETER
     +         (10,'conversion error of ASCII -> number ...')
C
         ELSE IF ((WORKY.EQ.'P').OR.(WORKY.EQ.'Q')) THEN
            LL = LL + 1
         ENDIF
250   CONTINUE
C      
C  extract 1-arg, 2-arg operations and save them (with offsets)
300   WORK1(1:) = ZLINE(1:)
      LL = (WORKLM+1) / 2
      DO 350, N=1,999
         CALL EXPRDC(WORK1,WORK,OPER(N),PPP(N))
         IF (PPP(N).LT.1) CALL STETER(1,ERROR3)
         IF (WORK(2:2).NE.' ') THEN
            WORK1 = WORK
         ELSE
            IF (N.GT.LL) CALL STETER
     +         (11,'too many operations in expression - split it up...')
            OPCNT = N 
            GOTO 1000
         ENDIF
350   CONTINUE
C 
1000  IF (ZDEFA(2:2).EQ.'X') THEN               !we must have a result frame
         CALL STFOPN(RESFRA,D_R8_FORMAT,0,F_IMA_TYPE,RIMNO,STAT)
         CALL STDRDI(RIMNO,'NAXIS',1,1,IAV,NAXIS,UNIT,NNULL,STAT)
         IF (NAXIS.LT.3) 
     +      CALL STETER(25,'invalid 2-dim result frame...')
         CALL STDRDI(RIMNO,'NPIX',1,3,IAV,NPIX,UNIT,NNULL,STAT)
         IF (RFELEM.GT.NPIX(3)) CALL STETER(24,ERROR5)
         SIZE = NPIX(1)*NPIX(2)
         RFELEM = (RFELEM-1)*SIZE+ 1               !now first elem
         NAXIS = 2
         CALL STDRDR(RIMNO,'LHCUTS',1,4,IAV,RCUTS,UNIT,NNULL,STAT)
         CUTS(3) = DBLE (RCUTS(3))
         CUTS(4) = DBLE (RCUTS(4))
      ENDIF
C      
C  if no frames involved ... do it right now + skip the loop
C      
      IF (K.LE.0) THEN
         CALL COMPUX(PNTR,OPER,PPP,OPCNT,APNTRS,ATOM,DATTA,PIDX,
     +               COUNT,NPIX(1),RPNTR,CUTS(3))
         CONST = CUTS(3)
         IF (ZDEFA(1:1).EQ.'N') THEN    !put resulting constant into OUTPUTD(1)
            CALL STKWRD('OUTPUTD',CONST,1,1,UNIT,STAT)
            WRITE(OUTPU,10000) CONST
            CALL STTPUT(OUTPU,STAT)
         ELSE
            IF (ZDEFA(2:2).NE.'X') THEN 
               CALL STIGET(RESFRA,D_R8_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +                     6,NAXIS,NPIX,START,STEP,IDENT,
     +                     CUNIT,RPNTR,RIMNO,STAT)
               CALL STDRDR(RIMNO,'LHCUTS',1,4,IAV,RCUTS,UNIT,NNULL,STAT)
               SIZE = 1
               DO 1100, N=1,NAXIS
                  SIZE = SIZE * NPIX(N)
1100           CONTINUE
               CALL DOPCF(CONST,MADRID(RPNTR),SIZE)    !copy constant to frame
               RCUTS(3) = REAL (CONST)
               RCUTS(4) = REAL (CONST)
            ELSE
               CALL STFCRE('WORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +                     SIZE,WIMNO,STAT)
               CALL STFMAP(WIMNO,F_X_MODE,1,SIZE,IAV,RPNTR,STAT)
               CALL DOPCF(CONST,MADRID(RPNTR),SIZE)    !copy constant to frame
               CALL STFPUT(RIMNO,RFELEM,SIZE,MADRID(RPNTR),STAT)
               IF (CONST.LT.CUTS(3)) THEN
                  CUTS(3) = CONST
               ELSEIF (CONST.GT.CUTS(4)) THEN
                  CUTS(4) = CONST
               ELSE
                  GOTO 9000
               ENDIF
               RCUTS(3) = REAL (CUTS(3))
               RCUTS(4) = REAL (CUTS(4))
            ENDIF
            CALL STDWRR(RIMNO,'LHCUTS',RCUTS,1,4,UNIT,STAT)
            CALL DSCUPT(RIMNO,RIMNO,' ',STAT)
         ENDIF
         GOTO 9000
      ENDIF
C      
C  at least one frame involved ...
C      
      FRACNT = K
      CALL ISSUBF(RESFRA,K)
      IF (K.GT.0) THEN
         UPDA = 1
      ELSE 
         UPDA = 0
      ENDIF
C 
      IF (DSCFRA(1:1).NE.' ') THEN
         DSCORG = 0
      ELSE
         DSCORG = 1
      ENDIF
C 
C  loop over all input frames to check, if we update result frame
C 
      DO 1250, N=1,FRACNT
         CALL STFINF(INFRAM(N),2,INFO,STAT)
         CALL STFOPN(INFRAM(N),D_R8_FORMAT,0,F_IMA_TYPE,IMNO(N),STAT)
C 
         IF (ZDEFA(2:2).NE.'X')THEN
            IF (UPDA.EQ.0) 
     +         CALL GENEQF(INFRAM(N),RESFRA,UPDA)   !input_frame = res_frame?
            IF (DSCORG.EQ.0)  THEN
               CALL GENEQF(INFRAM(N),DSCFRA,LL)
               IF (LL.NE.0) DSCORG = N
            ENDIF
         ENDIF
1250  CONTINUE
C 
      IF (ZDEFA(2:2).EQ.'X')THEN
         UPDA = 1
         DSCORG = 0
      ELSE IF (DSCORG.EQ.0) THEN
         DSCORG = 1                                 !if not found take 1. frame
      ENDIF
C 
C  if no frame is equal to result frame, we create it
C
      IF (UPDA.EQ.0) THEN                !get standard descr. from DSCORG frame
         CALL STDRDI(IMNO(DSCORG),
     +               'NAXIS',1,1,IAV,NAXIS,UNIT,NNULL,STAT)
         CALL STDRDI(IMNO(DSCORG),
     +               'NPIX',1,NAXIS,IAV,NPIX,UNIT,NNULL,STAT)
         CALL STDRDD(IMNO(DSCORG),
     +               'START',1,NAXIS,IAV,START,UNIT,NNULL,STAT)
         CALL STDRDD(IMNO(DSCORG),
     +               'STEP',1,NAXIS,IAV,STEP,UNIT,NNULL,STAT)
         CALL STDRDC(IMNO(DSCORG),
     +               'IDENT',1,1,72,IAV,IDENT,UNIT,NNULL,STAT)
         M = (NAXIS+1) * 16
         CALL STDRDC(IMNO(DSCORG),
     +               'CUNIT',1,1,M,IAV,CUNIT,UNIT,NNULL,STAT)
C      
         SIZE = 1
         DO 1300, N=1,NAXIS
            SIZE = SIZE*NPIX(N)
1300     CONTINUE
         CALL STFCRE(RESFRA,DATFMT,F_O_MODE,F_IMA_TYPE,
     +               SIZE,RIMNO,STAT)
         CALL STDWRI(RIMNO,'NAXIS',NAXIS,1,1,UNIT,STAT)
         CALL STDWRI(RIMNO,'NPIX',NPIX,1,NAXIS,UNIT,STAT)
         CALL STDWRD(RIMNO,'START',START,1,NAXIS,UNIT,STAT)
         CALL STDWRD(RIMNO,'STEP',STEP,1,NAXIS,UNIT,STAT)
         CALL STDWRC(RIMNO,'IDENT',1,IDENT,1,72,UNIT,STAT)
         CALL STDWRC(RIMNO,'CUNIT',1,CUNIT,1,M,UNIT,STAT)
C 
C  else we just open the frame
C 
      ELSE
         IF (ZDEFA(2:2).NE.'X') THEN
            CALL STFOPN(RESFRA,DATFMT,0,F_IMA_TYPE,RIMNO,STAT)
            CALL STDRDI(RIMNO,'NAXIS',1,1,IAV,NAXIS,UNIT,NNULL,STAT)
            CALL STDRDI(RIMNO,'NPIX',1,NAXIS,IAV,NPIX,UNIT,NNULL,STAT)
            SIZE = 1
            DO 1330, N=1,NAXIS
               SIZE = SIZE*NPIX(N)
1330        CONTINUE
         ENDIF
         CALL STFINF(RESFRA,2,INFO,STAT)
         IF (INFO(2).NE.D_R8_FORMAT) THEN
            CUNIT(1:) =
     +      'Warning: Result frame is not real - '//
     +      'you may get truncation errors... '
            CALL STTPUT(CUNIT,STAT)
         ENDIF
         CALL STDRDD(RIMNO,'START',1,NAXIS,IAV,START,UNIT,NNULL,STAT)
         CALL STDRDD(RIMNO,'STEP',1,NAXIS,IAV,STEP,UNIT,NNULL,STAT)
         DSCORG = 0
      ENDIF
C 
C  here we determine, if we have to use windows or take the complete frames
C  also COMPUTE/PIX command or if just one frame ...
C      
      IF (ZDEFA(2:2).EQ.'X') THEN
         DO 1600, N=1,FRACNT
            CALL STDRDI(IMNO(N),
     +                  'NAXIS',1,1,IAV,NAXISA,UNIT,NNULL,STAT)
            CALL STDRDI(IMNO(N),
     +                  'NPIX',1,NAXISA,IAV,NPIXA,UNIT,NNULL,STAT)
            IF (FELEM(N).GT.1) THEN
               IF ((NAXISA.LT.3) .OR. (NPIXA(3).LT.FELEM(N)))
     +            CALL STETER(27,ERROR5)
               FELEM(N) = (FELEM(N)-1)*SIZE + 1
            ELSE IF (NAXISA.LT.2)THEN
               CALL STETER(26,'1-dim input frame not supported...')
            ENDIF
            DO 1500, K=1,NAXIS
               IF (NPIX(K).NE.NPIXA(K)) 
     +            CALL STETER(26,'non-matching NPIX...')
1500        CONTINUE
1600     CONTINUE
C 
      ELSE IF (ZDEFA(2:2).EQ.'P') THEN
         DO 1800, N=1,FRACNT
            IF (N.NE.DSCORG) THEN
               CALL STDRDI(IMNO(N),
     +                     'NAXIS',1,1,IAV,NAXISA,UNIT,NNULL,STAT)
               IF (NAXISA.NE.NAXIS)
     +            CALL STETER(26,'non-matching NAXIS...')
               CALL STDRDI(IMNO(N),
     +                     'NPIX',1,NAXISA,IAV,NPIXA,UNIT,NNULL,STAT)
               DO 1700, K=1,NAXIS
                  IF (NPIX(K).NE.NPIXA(K)) 
     +               CALL STETER(26,'non-matching NPIX...')
1700           CONTINUE
            ENDIF
1800     CONTINUE
C      
      ELSE
         DO 2000, N=1,FRACNT
            IF (N.NE.DSCORG) THEN
               CALL STDRDI(IMNO(N),
     +                     'NPIX',1,NAXIS,IAV,NPIXA,UNIT,NNULL,STAT)
               CALL STDRDD(IMNO(N),
     +                     'START',1,NAXIS,IAV,STARTA,UNIT,NNULL,STAT)
               CALL STDRDD(IMNO(N),
     +                     'STEP',1,NAXIS,IAV,STEPA,UNIT,NNULL,STAT)
               DO 1900, K=1,NAXIS
                  IF (NPIX(K).NE.NPIXA(K)) GOTO 3000
                  EPS = 0.001 * ABS(STEP(K))                  !0.1 % of step
                  DIF = ABS(STEP(K) - STEPA(K))
                  IF (DIF.GT.EPS) GOTO 3000
                  EPS = 0.01 * ABS(STEP(K))                !1.0 % of step
                  DIF = ABS(START(K) - STARTA(K))
                  IF (DIF.GT.EPS) GOTO 3000
1900           CONTINUE
            ENDIF
2000     CONTINUE
         GOTO 5000
C      
C  frames are not of same size or unaligned
C  so we only work on overlapping area (done in COMPUY)
C  but open DESCR originator frame again, so we can use DSCUPT after
C  
3000     CALL COMPUY(DATFMT,ZLINE,ATOM,APNTRS,UPDA,RESFRA,CUTS(3))
         CALL STFOPN(RESFRA,DATFMT,0,F_IMA_TYPE,RIMNO,STAT)
         IF (DSCORG.NE.0) THEN
            CALL STFOPN(INFRAM(DSCORG),
     +                  DATFMT,0,F_IMA_TYPE,IMNO(DSCORG),STAT)
            CALL DSCUPT(IMNO(DSCORG),RIMNO,' ',STAT)
         ELSE
            CALL DSCUPT(RIMNO,RIMNO,' ',STAT)
         ENDIF
         RCUTS(3) = REAL(CUTS(3))
         RCUTS(4) = REAL(CUTS(4))
         CALL STDWRR(RIMNO,'LHCUTS',RCUTS,1,4,UNIT,STAT)
         GOTO 9000
      ENDIF
C      
C  here we can work in pixel space
C      
5000  CALL STKRDI('MONITPAR',20,1,LL,MAPSIZ,UNIT,NNULL,STAT)
      IF (MAPSIZ.LT.100) MAPSIZ = 100               !default chunk size
      MAPSIZ = MAPSIZ*MAPSIZ
      IF (FRACNT.GT.3) MAPSIZ = MAPSIZ/2
      IF (MAPSIZ.GT.SIZE) THEN
         MAPSIZ = SIZE
         TOTAL = 1
      ELSE
         TOTAL = SIZE / MAPSIZ
         IF ( (TOTAL*MAPSIZ) .LT. SIZE) TOTAL = TOTAL + 1
      ENDIF
C 
C  map virtual memory for input files, 4 working buffers + result frame
      K = (FRACNT + 5) * MAPSIZ
      CALL STFCRE('WORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            K,WIMNO,STAT)
      CALL STFMAP(WIMNO,F_X_MODE,1,K,IAV,PNTR(1),STAT)
      IF (FRACNT.GT.1) THEN
         DO 5200, N=2,FRACNT                 !no more pointers,
            PNTR(N) = PNTR(N-1) + 2*MAPSIZ     !but indices of MADRID ...
5200     CONTINUE
      ENDIF
C 
      RPNTR = PNTR(FRACNT) + 2*MAPSIZ         !double precision...
C 
C  get pointers for 4 working chunks
      DO 5250, N=1,4
         PNTR(N+23) = RPNTR + (2*N*MAPSIZ)
5250  CONTINUE
C      
C  loop over all chunks
C
      DO 5800, N=1,TOTAL
         DO 5600, K=1,FRACNT
            CALL STFGET(IMNO(K),FELEM(K),MAPSIZ,TRSIZE,
     +                  MADRID(PNTR(K)),STAT)
            FELEM(K) = FELEM(K) + TRSIZE
5600     CONTINUE
C      
         CALL COMPUX(PNTR,OPER,PPP,OPCNT,APNTRS,ATOM,DATTA,PIDX,
     +               COUNT,TRSIZE,RPNTR,CUTS(3))
C      
         CALL STFPUT(RIMNO,RFELEM,TRSIZE,MADRID(RPNTR),STAT)
         RFELEM = RFELEM + TRSIZE
5800  CONTINUE
C  
C  ***
C  we're done - set dynamic range of result frame
C  ***
C      
C  take care of descriptors
      IF (DSCORG.NE.0) THEN
         CALL DSCUPT(IMNO(DSCORG),RIMNO,' ',STAT)
      ELSE IF (UPDA.EQ.1) THEN                       !just update history
         CALL DSCUPT(RIMNO,RIMNO,' ',STAT)
      ENDIF
      RCUTS(3) = REAL(CUTS(3))
      RCUTS(4) = REAL(CUTS(4))
      CALL STDWRR(RIMNO,'LHCUTS',RCUTS,1,4,UNIT,STAT)
C      
C  That's it folks...
9000  RETURN
C
C  formats...
10000 FORMAT(G24.15)
      END

      SUBROUTINE COMPUX(PNTR,OPER,PPP,OPCNT,APNTRS,ATOM,DATTA,PIDX,
     +                     COUNT,MAPSIZ,RPNTR,CUTS)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine COMPUX			951011
C  K. Banse				ESO - Garching
C
C.KEYWORDS
C  arithmetic operations
C
C.PURPOSE
C  evaluate the arithmetic expression specified in the main module
C  for the currently mapped section of the complete data space
C
C.ALGORITHM
C  straight forward but rather tricky - we have to juggle pointers like hell...
C
C.INPUT/OUTPUT
C  call as COMPUX(PNTR,OPER,PPP,OPCNT,APNTRS,ATOM,DATTA,PIDX,
C               COUNT,MAPSIZ,RPNTR,CUTS)
C  input par:
C  PNTR:	I*8 array	data pointers
C  OPER:	char. array	holds the individual operations
C  				which make up the complete expression
C  PPP:		I*4 array	points to the first atom to be used 
C  				with OPER(n)
C  OPCNT:	I*4		size of OPER + PPP
C  APNTRS:	I*4 array	initial array of pointers to ATOM
C  ATOM:	char. array 	holds names of frames, names of functions
C  				and constants in ASCII
C  				only needed for function names...
C  DATTA:	R*8 array	holds constants already as reals
C  PIDX:	I*4 array	array of pointers to frames
C  COUNT:	I*4		total no. of operands in original expression
C  MAPSIZ:	I*4		size of mapped chunks
C  RPNTR:	I*8		pointer to result frame
C  
C  output par:
C  CUTS:	R*8 array	minmax values for result frame
C      
C.VERSIONS
C  see SCCS
C      
C-------------------------------------------------------------------------
C      
      IMPLICIT NONE
C
      CHARACTER      ATOM(23)*80,OPER(23)*5
      CHARACTER      BTOM(23)*80
      CHARACTER      OPERA*5
C      
      DOUBLE PRECISION  DATTA(*)
      DOUBLE PRECISION  CUTS(2),CNSTAB(2),CONST
      DOUBLE PRECISION  BDATA(23),FMIN,FMAX
C      
      INTEGER      PPP(*),OPCNT,APNTRS(*),PIDX(*),COUNT
      INTEGER      BPNTRS(48),BPIDX(23)
      INTEGER*8    RPNTR,OUTPTR,KPNTR,LPNTR,PNTR(27)
      INTEGER      MAPSIZ
      INTEGER      FRACNT
      INTEGER      N,IDUMMY
      INTEGER      NBRA,P1,P2,P3,PP
      INTEGER      OPINDX,NOUT
      INTEGER      MADRID(1)
C      
      COMMON      /IPOOL/  FRACNT
C 
      COMMON      /VMR/ MADRID
C      
C  first initialize pointers for ATOM
      OPINDX = 0
      DO 100 N=1,COUNT
         BTOM(N) = ATOM(N)
         BDATA(N) = DATTA(N)
         BPIDX(N) = PIDX(N)
100   CONTINUE
      FMIN = CUTS(1)
      FMAX = CUTS(2)
C 
      DO 150, N=1,48
         BPNTRS(N) = APNTRS(N)
150   CONTINUE
      IDUMMY = 0                        !init DUMMY counter
C      
C  extract basic OPERAs
1000  OPINDX = OPINDX + 1
      OPERA = OPER(OPINDX)
      PP = PPP(OPINDX)
C      
C  extract operands
      P1 = BPNTRS(PP)
      P2 = BPNTRS(PP+1)
      P3 = BPNTRS(PP+2)
C      
C  find out what kind of OPERA to do
      IF (INDEX(OPERA,'C').GT.0) THEN
         IF (INDEX(OPERA,'F').GT.0) GOTO 1600            !also F involved...?
      ELSE
         IF (OPERA(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 (OPERA(1:2).EQ.'CC') THEN
         CALL DOPCC(OPERA,BDATA(P1),BDATA(P2),CONST)
      ELSE
         IF (OPERA(1:1).EQ.'P') THEN                  !1-arg functions
            CALL DF1C(BTOM(P1)(1:5),BDATA(P2),CONST)
         ELSE                                          !2-arg functions
            CNSTAB(1) = BDATA(P2)
            CNSTAB(2) = BDATA(P3)
            CALL DF2CC(BTOM(P1)(1:5),CNSTAB,CONST)
         ENDIF
      ENDIF
C
C  put resulting constant back into relevant BTOM + goto loopend
      IF (OPCNT.NE.OPINDX) THEN
         BDATA(P1) = CONST
         GOTO 5000
      ELSE                        !we're done already...
         CUTS(1) = CONST
         GOTO 8000
      ENDIF
C
C  ***
C  one operand is a file
C  ***
C      
1600  IF (OPERA(1:2).EQ.'FC') THEN
         KPNTR = PNTR(BPIDX(P1))
         CONST = BDATA(P2)
         NBRA = 1                                    !no function
      ELSE IF (OPERA(1:2).EQ.'CF') THEN
         KPNTR = PNTR(BPIDX(P2))
         CONST = BDATA(P1)
         NBRA = 1                                    !no function
      ELSE IF (OPERA(1:2).EQ.'PF') THEN
         KPNTR = PNTR(BPIDX(P2))
         NBRA = 2                                    !1-arg function
      ELSE
         NBRA = 3                                    !2-arg function
         IF (OPERA(1:2).NE.'QC') THEN
            KPNTR = PNTR(BPIDX(P2))
            CONST = BDATA(P3)
         ELSE
            KPNTR = PNTR(BPIDX(P3))
            CONST = BDATA(P2)
         ENDIF
      ENDIF
C      
C determine where result should go to
C      
      IF (OPCNT.EQ.OPINDX) THEN            !last OPERA ?
         OUTPTR = RPNTR
      ELSE
         DO 1700, N=24,27                  !24,25,26,27 are the dummy pointers
            IF (KPNTR.EQ.PNTR(N)) THEN     !operand is one of the dummy frames
               OUTPTR = PNTR(N)
               NOUT = N
               GOTO 1900
            ENDIF
1700     CONTINUE
C                                    !we have to use up one of the dummy frames
         IF (IDUMMY.LE.3) THEN
            IDUMMY = IDUMMY + 1
            NOUT = 23 + IDUMMY            !keep no. of dummy area used 
            OUTPTR = PNTR(NOUT)            !for output
         ELSE
            CALL STETER(99,'error in internal logic ...!')
         ENDIF
      ENDIF
C
C  now do the actual OPERA
1900  IF (NBRA.EQ.1) THEN
         CALL DOPFC(OPERA,MADRID(KPNTR),CONST,
     +                 MADRID(OUTPTR),MAPSIZ)
      ELSE IF (NBRA.EQ.2) THEN
         CALL DF1F(BTOM(P1)(1:5),MADRID(KPNTR),
     +                  MADRID(OUTPTR),MAPSIZ)
      ELSE
         CALL DF2FC(BTOM(P1)(1:5),MADRID(KPNTR),CONST,
     +                   MADRID(OUTPTR),MAPSIZ)
      ENDIF
C      
      IF (OPCNT.NE.OPINDX) THEN
         GOTO 4400                              !goto loopend
      ELSE
         GOTO 7700
      ENDIF
C
C  ***
C  both operands are files 
C  ***
C      
2000  IF (OPERA(1:1).NE.'Q') THEN
         KPNTR = PNTR(BPIDX(P1))
         LPNTR = PNTR(BPIDX(P2))
      ELSE
         KPNTR = PNTR(BPIDX(P2))
         LPNTR = PNTR(BPIDX(P3))
      ENDIF
C      
C determine where result should go to
C      
      IF (OPCNT.EQ.OPINDX) THEN            !last OPERA ?
         OUTPTR = RPNTR
      ELSE
         DO 2100, N=24,27                  !24,25,26,27 are the dummy pointers
            IF ( (KPNTR.EQ.PNTR(N)) .OR.   !operand is one of the dummy frames
     +           (LPNTR.EQ.PNTR(N)) ) THEN      
               OUTPTR = PNTR(N)
               NOUT = N
               GOTO 4000
            ENDIF
2100     CONTINUE
C                              !we have to use up one of the dummy frames
         IF (IDUMMY.LE.3) THEN
            IDUMMY = IDUMMY + 1
            NOUT = 23 + IDUMMY            !keep no. of dummy area used 
            OUTPTR = PNTR(NOUT)            !for output
         ELSE
            CALL STETER(99,'error in internal logic ...!')
         ENDIF
      ENDIF
C
C  now do the actual OPERA
4000  IF (OPERA(1:1).NE.'Q') THEN
         CALL DOPFF(OPERA,MADRID(KPNTR),MADRID(LPNTR),
     +                   MADRID(OUTPTR),MAPSIZ)
      ELSE
         CALL DF2FF(BTOM(P1)(1:5),MADRID(KPNTR),MADRID(LPNTR),
     +                   MADRID(OUTPTR),MAPSIZ)
      ENDIF
C      
      IF (OPCNT.EQ.OPINDX) GOTO 7700
C
C  put resulting frame pointer back into relevant BPIDX, if we are not finished yet
4400  BPIDX(P1) = NOUT
C      
C  loopend for all basic OPERAs
5000  IF (OPERA(1:1).NE.'Q') THEN
         NBRA = 2
      ELSE
         NBRA = 3
      ENDIF
      DO 5100, N=PP+1,45                        !update pointers for BTOM
         BPNTRS(N) = BPNTRS(N+NBRA)
5100  CONTINUE
      GOTO 1000                        !get next OPERA
C      
C  finally get min + max
7700  CALL DMYMX(MADRID(RPNTR),MAPSIZ,CUTS)
      IF (CUTS(1).GT.FMIN) CUTS(1) = FMIN
      IF (CUTS(2).LT.FMAX) CUTS(2) = FMAX
C      
C  that's it folks...
8000  RETURN
C      
      END