C @(#)arthmy.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 ARTHMY(PLADIR,ZDEFA,ZRESFR,ZLINE)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine ARTHMY			version 1.00	950803
C  K. Banse            		      	ESO - Garching
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  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 ARTHMY(PLADIR,ZDEFA,ZRESFR,ZLINE)
C  
C  input par:
C  PLADIR:	I*4		plane direction 
C				X_Z_PLANE, X_Y_PLANE or Z_Y_PLANE
C  ZDEFA:	char.		2 char. default string
C				(1) Y or N for result frame or not
C  ZRESFR:	char.		name of result frame
C  ZLINE:	char.		char.string with arithmetic expression
C	
C.VERSIONS
C  see SCCS
C	
C 000621
C-------------------------------------------------------------------------
C
      IMPLICIT NONE
C
C  important parameter: WORKLM = 45
C
      CHARACTER*(*)  ZDEFA,ZRESFR,ZLINE
      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*40
C      
      REAL         CONST,CUTS(4)
      REAL         DATTA(23),RDUM(1)                 !size = MAXFRAM
C 
      DOUBLE PRECISION STEP(3),START(3),DDUM(1)
C      
      INTEGER      LATOM(23),WORKLM                  !size = MAXFRAM
      INTEGER      APNTRS(48)                        !size = WORKLM + 3
      INTEGER*8    PNTR(27)                          !size = MAXFRAM + 4
      INTEGER*8    XPNTR,RPNTR
      INTEGER      IMNO(23),IPLANE(23)               !size = MAXFRAM
      INTEGER      DIM3(23),PIDX(23),PPP(23)         !size = MAXFRAM
      INTEGER      PLADIR,OPLANE,PLINDX
      INTEGER      FRACNT,NPIX(3),NPIXA(3),YPIX(2)
      INTEGER      COUNT,IAV,K,LL,DIM1ON
      INTEGER      N,NAXIS,NAXISA,TOTAL
      INTEGER      SIZE,INFO(5),DATFMT
      INTEGER      MAPSIZ,STAT
      INTEGER      NNULL,UNIT(1),INPU(1),OPCNT
      INTEGER      RIMNO,WIMNO,MADRID(1)
      INTEGER      UPDA,DSCORG,XIMNO
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C      
      COMMON      /IPOOL/  FRACNT
      COMMON      /VMR/ MADRID
C 
      DATA      WORKLM      /45/
      DATA      CUNIT       /' '/,  IDENT       /' '/
      DATA      ERROR3
     +          /'wrong syntax in arithmetic expression ...'/
      DATA      ERROR4      /'too many operands...'/
C      
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C      
C  init variables
      TOTAL = 0
      MAPSIZ = 0
      DIM1ON = 0
      CUTS(1) = 0.
      CUTS(2) = 0.
      CUTS(3) = + 999999.9
      CUTS(4) = - 999999.9
C 
      DO 130, N=1,23
         DATTA(N) = 0.0
130   CONTINUE
      DO 140, N=1,48
         APNTRS(N) = -1
140   CONTINUE
      DO 150, N=1,3
         NPIX(N) = 1
         START(N) = 0.D0
         STEP(N) = 1.D0
150   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_R4_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 because 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
      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
            CALL CLNFRA(ATOM(LL),INFRAM(K),0)
            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 UPCAS(ATOM(LL),ATOM(LL))
            STAT = INDEX(ATOM(LL),'D')
            IF (STAT.GT.1) THEN
               DATFMT = D_R8_FORMAT     !explicit double constant...
               CALL GENCNV(ATOM(LL),4,1,INPU,RDUM,DDUM,STAT)
               DATTA(LL) = SNGL(DDUM(1))
            ELSE
               CALL GENCNV(ATOM(LL),2,1,INPU,RDUM,DDUM,STAT)
               DATTA(LL) = RDUM(1)
            ENDIF
            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      
C  if no frames involved ... do it right now + skip the loop
C      
1000  IF (K.LE.0) THEN
         CALL COMPUA(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 OUTPUTR(1)
            CALL STKWRR('OUTPUTR',CONST,1,1,UNIT,STAT)
            WRITE(OUTPU,10000) CONST
            CALL STTPUT(OUTPU,STAT)
         ELSE
            CALL STIGET(RESFRA,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +                  3,NAXIS,NPIX,START,STEP,IDENT,
     +                  CUNIT,RPNTR,RIMNO,STAT)
            SIZE = NPIX(1) * NPIX(2) * NPIX(3)
            CALL CONFIL(CONST,MADRID(RPNTR),SIZE)       !copy constant to frame
            CUTS(4) = CONST
         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
      IF (DSCFRA(1:1).NE.' ') THEN
         DSCORG = 0
      ELSE
         DSCORG = 1
      ENDIF
C 
      DO 1200, N=1,FRACNT
         CALL STFINF(INFRAM(N),2,INFO,STAT)
         IF (INFO(2).EQ.D_R8_FORMAT) DATFMT = D_R8_FORMAT
         CALL STFOPN(INFRAM(N),D_R4_FORMAT,0,F_IMA_TYPE,IMNO(N),STAT)
         IF (UPDA.EQ.0) THEN
            CALL GENEQF(INFRAM(N),RESFRA,UPDA)       !test, if frames are equal
            IF (UPDA.NE.0) UPDA = N                  !Yes. Save the index
         ENDIF
         IF (DSCORG.EQ.0)  THEN
            CALL GENEQF(INFRAM(N),DSCFRA,LL)
            IF (LL.NE.0) DSCORG = N
         ENDIF
1200  CONTINUE
C 
      IF (PLADIR.EQ.X_Y_PLANE) THEN
         PLINDX = 3                              !index of `looping' plane 
      ELSE IF (PLADIR.EQ.X_Z_PLANE) THEN
         PLINDX = 2
      ELSE
         PLINDX = 1
      ENDIF
C 
C  we use R*4 data currently
      IF (DATFMT.EQ.D_R8_FORMAT) THEN
         CUNIT(1:) = 
     +   'Warning: double operand found - result will be single! '
         CALL STTPUT(CUNIT,STAT)
         DATFMT = D_R4_FORMAT
      ENDIF
C 
C  work on first file operand 
      CALL STDRDI(IMNO(1),'NAXIS',1,1,IAV,INPU,UNIT,NNULL,STAT)
      NAXIS = INPU(1)
      CALL STDRDI(IMNO(1),'NPIX',1,NAXIS,IAV,NPIX,UNIT,NNULL,STAT)
      IF (NAXIS.LT.2) THEN
         DIM3(1) = -1
         DIM1ON = 1                        !so we know later on...
      ELSE IF (NAXIS.EQ.2) THEN
         IPLANE(1) = X_Y_PLANE
         TOTAL = 1
         DIM3(1) = 1
         MAPSIZ = NPIX(1) * NPIX(2)
      ELSE
         IPLANE(1) = PLADIR
         TOTAL = NPIX(PLINDX)
         DIM3(1) = TOTAL
         SIZE = NPIX(1) * NPIX(2) * NPIX(3)
         MAPSIZ = SIZE / NPIX(PLINDX)
         CALL STDRDD(IMNO(1),'START',1,NAXIS,IAV,START,UNIT,NNULL,STAT)
         CALL STDRDD(IMNO(1),'STEP',1,NAXIS,IAV,STEP,UNIT,NNULL,STAT)
         CALL STDRDC(IMNO(1),'IDENT',1,1,72,IAV,IDENT,UNIT,NNULL,STAT)
         CALL STDRDC(IMNO(1),'CUNIT',1,1,64,IAV,CUNIT,UNIT,NNULL,STAT)
         IF (DSCORG.EQ.0) DSCORG = 1
      ENDIF
C 
C  loop over all file operands
      IF (FRACNT.GT.1) THEN
         DO 1500, N=2,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 (NAXISA.LT.2) THEN
               DIM3(N) = -1
               DIM1ON = 1                        !so we know later on...
               GOTO 1500                         !avoid size check
            ELSE IF (NAXISA.EQ.2) THEN
               IPLANE(N) = X_Y_PLANE
               DIM3(N) = 1
               K = NPIXA(1) * NPIXA(2)
            ELSE
               IPLANE(N) = PLADIR
               DIM3(N) = NPIXA(PLINDX)
               IF (TOTAL.LE.1) THEN
                  TOTAL = DIM3(N)
                  NPIX(1) = NPIXA(1)
                  NPIX(2) = NPIXA(2)
                  NPIX(3) = NPIXA(3)
                  SIZE = NPIX(1) * NPIX(2) * NPIX(3)
                  CALL STDRDD(IMNO(N),'STEP',1,NAXIS,IAV,START,UNIT,
     +                        NNULL,STAT)
                  CALL STDRDD(IMNO(N),'STEP',1,NAXIS,IAV,STEP,UNIT,
     +                        NNULL,STAT)
                  CALL STDRDC(IMNO(N),'IDENT',1,1,72,IAV,IDENT,UNIT,
     +                        NNULL,STAT)
                  CALL STDRDC(IMNO(N),'CUNIT',1,1,64,IAV,CUNIT,UNIT,
     +                        NNULL,STAT)
                  IF (DSCORG.EQ.0) DSCORG = N
               ELSE IF (DIM3(N).NE.TOTAL) THEN
                  CALL STETER (67,
     +            'No. of planes must be the same for all frames')
               ENDIF
               K = (NPIXA(1) * NPIXA(2) * NPIXA(3)) / NPIXA(PLINDX)
            ENDIF
            IF (MAPSIZ.EQ.0) THEN
               MAPSIZ = K
            ELSE
               IF (K.NE.MAPSIZ) THEN
                  CALL STETER(66,
     +            'All relevant planes must have equal size')
               ENDIF
            ENDIF
1500     CONTINUE
      ENDIF
C      
      YPIX(1) = NPIX(1)                  !maybe we need it 
      YPIX(2) = NPIX(2)
      IF (TOTAL.EQ.0) THEN               !only 1-dim frames involved
         CALL STETER(68,
     +   'Only 1-dim files involved - use COMPUTE/IMAGE instead.')
      ENDIF
      IF (TOTAL.EQ.1) THEN               !only 2-dim frames involved
         SIZE = NPIX(1) * NPIX(2)
         OPLANE = X_Y_PLANE
         NAXIS = 2
      ELSE
         OPLANE = PLADIR
         NAXIS = 3
         IF (DIM1ON.EQ.1) THEN           !we have 1-dime operands
            IF (OPLANE.EQ.X_Z_PLANE) THEN
               YPIX(2) = NPIX(3)
            ELSE IF (OPLANE.EQ.Z_Y_PLANE) THEN
               YPIX(1) = NPIX(3)
               YPIX(2) = NPIX(2)
            ENDIF                        !YPIX for X_Y_PLANE already set
         ENDIF 
      ENDIF
C 
C  for 1-dim operands we need one more dimension check
      IF (DIM1ON.EQ.1) THEN
         DO 2000, N=1,FRACNT
            IF (DIM3(N).EQ.-1) THEN
               CALL STDRDI(IMNO(N),'NPIX',1,1,IAV,NPIXA,UNIT,
     +                     NNULL,STAT)
               IF (NPIXA(1).NE.YPIX(1)) THEN
                  CALL STETER(69,
     +       '1-dim frames must have same size as first plane dim.')
               ENDIF
            ENDIF
2000     CONTINUE
C 
         CALL STFCRE('Work1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               YPIX(1),XIMNO,STAT)
         CALL STFMAP(XIMNO,F_X_MODE,1,YPIX(1),IAV,XPNTR,STAT)
      ENDIF
C      
C  if no frame is equal to result frame, we create it
      IF (UPDA.EQ.0) THEN
         CALL STFCRE(RESFRA,DATFMT,F_O_MODE,F_IMA_TYPE,
     +               SIZE,RIMNO,STAT)
         INPU(1) = NAXIS
         CALL STDWRI(RIMNO,'NAXIS',INPU,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)
         N = (NAXIS+1)*16
         CALL STDWRC(RIMNO,'CUNIT',1,CUNIT,1,N,UNIT,STAT)
      ELSE
         IF (IPLANE(UPDA).NE.OPLANE) THEN
            CALL STETER(69,
     +  'We cannot update result frame with different plane directions')
         ENDIF
         CALL STFOPN(RESFRA,DATFMT,0,F_IMA_TYPE,RIMNO,STAT)
         CALL STFINF(RESFRA,2,INFO,STAT)
         IF ( (INFO(2).NE.D_R8_FORMAT) .AND.
     +        (INFO(2).NE.D_R4_FORMAT) ) THEN
            CUNIT(1:) = 
     +      'Warning: Result frame is not real - '//
     +      'you may get truncation errors... '
            CALL STTPUT(CUNIT,STAT)
         ENDIF
      ENDIF
C 
C  map virtual memory for input files + result frame
      K = (FRACNT + 5) * MAPSIZ
      CALL STFCRE('Work2',D_R4_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) + MAPSIZ      !but indices of MADRID ...
5200     CONTINUE
      ENDIF
C 
      RPNTR = PNTR(FRACNT) + MAPSIZ            !result pointer
C 
C  get pointers for 4 working chunks
      DO 5250, N=1,4
         PNTR(N+23) = RPNTR + (N*MAPSIZ)
5250  CONTINUE
C      
C  loop over all planes
C
      DO 5800, N=1,TOTAL
         DO 5600, K=1,FRACNT
            IF (DIM3(K).GE.N) THEN
               CALL STPGET(IMNO(K),IPLANE(K),N,MADRID(PNTR(K)),STAT)
            ELSE IF (DIM3(K).EQ.-1) THEN           !build 2-dim frame
               CALL STFGET(IMNO(K),1,YPIX(1),IAV,MADRID(XPNTR),STAT)
               CALL GROWIT('L',YPIX,MADRID(XPNTR),YPIX,
     +                     MADRID(PNTR(K)),1)
               DIM3(K) = 1
            ENDIF
5600     CONTINUE
C      
         CALL COMPUA(PNTR,OPER,PPP,OPCNT,APNTRS,ATOM,DATTA,PIDX,
     +               COUNT,MAPSIZ,RPNTR,CUTS(3))
C      
         CALL STPPUT(RIMNO,OPLANE,N,MADRID(RPNTR),STAT)
5800  CONTINUE
C  
C  ***
C  we're done - set dynamic range of result frame
C  ***
C      
C  take care of descriptors
8000  IF (FRACNT.GE.1) THEN
         CALL DSCUPT(IMNO(DSCORG),RIMNO,' ',STAT)
      ELSE
         CALL DSCUPT(RIMNO,RIMNO,' ',STAT)
      ENDIF
      CALL STDWRR(RIMNO,'LHCUTS',CUTS,1,4,UNIT,STAT)
C      
C  That's it folks...
9000  RETURN
C
C  formats...
10000 FORMAT(G15.7)
      END