* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM XVIEW_MAIN
C----------------------------------------------------------------------
C X-Window viewer for GILDAS images
C     F.Badia - F.Viallefond - S.Guilloteau
C----------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:errcod.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'mview.inc'
*
      CHARACTER*80 NAMEX
      CHARACTER*256 UIDDIR, HELPDIR
      INTEGER*4 POINTER
      INTEGER*4 BLC(4), TRC(4), PIDS(3)
      INTEGER I,FORM,N, LENC
      REAL BLANK(2)
      LOGICAL ERROR
      EXTERNAL USER_COORD
*
      DATA BLC /4*0/, TRC /4*0/
*
      CALL GILDAS_OPEN
      CALL GILDAS_CHAR('X_NAME$',NAMEX)
      CALL GILDAS_INTE('BLC$',BLC,4)
      CALL GILDAS_INTE('TRC$',TRC,4)
      CALL GILDAS_INTE('PROCESS_IDS$',PIDS,3)
      CALL GILDAS_CLOSE
*
      N = LENC(NAMEX)
      CALL SIC_PARSEF(NAMEX(1:N),X_FILE,' ','.gdf')
      CALL GDF_GEIS (X_ISLO,ERROR)
      IF (.NOT.ERROR) CALL
     $GDF_REIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-XVIEW,  Cannot read input file'
         CALL SYSEXI(FATALE)
      ENDIF
      CALL GDF_READX (X_ISLO,ERROR)
*
* Define projection and axis types.
      IF (X_PROJ.NE.0 .AND. X_PTYP.NE.0) THEN
         IF (X_XAXI.EQ.0.AND.X_YAXI.EQ.0) THEN
            X_XAXI = 1
            X_YAXI = 2
         ENDIF
         CALL GR8_PROJEC(X_PTYP,X_A0)
      ELSE
         X_XAXI = 1
         X_YAXI = 2
         CALL GR8_PROJEC(0,X_A0)
      ENDIF
      IF (X_FAXI.EQ.0) THEN
         DO I=1,3
            IF (I.NE.X_XAXI .AND. I.NE.X_YAXI) X_FAXI = I
         ENDDO
      ENDIF
*
* Blanking
      IF (X_BLAN.EQ.0) THEN
         X_BVAL = 0.0
         X_EVAL = -1.0
      ENDIF
*
* Define and load subset
      DO I=1,4
         BLC(I) = MAX(1,MIN(BLC(I),X_DIM(I)))
         IF (TRC(I).EQ.0) THEN
            TRC(I) = X_DIM(I)
         ELSE
            TRC(I) = MIN(X_DIM(I),TRC(I))
         ENDIF
         DIM(I) = TRC(I)-BLC(I)+1
      ENDDO
      OX = BLC(1)
      OY = BLC(2)
      OZ = BLC(3)
      FORM = FMT_R4
      CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,FORM,ERROR)
      IPX = POINTER(X_ADDR,MEMORY)
*
* Give to MVIEW
      N = LENC(X_FILE)+1
      X_FILE(N:N) = CHAR(0)
      UIDDIR = 'GILDAS_UID:'//CHAR(0)
      HELPDIR = 'GILDAS_RUN:'//CHAR(0)
      CALL SIC_GTLGTR(UIDDIR)
      CALL SIC_GTLGTR(HELPDIR)
      IF (UIDDIR.EQ.'GILDAS_UID:') UIDDIR = ' '
      N = LENC(UIDDIR)+1
      UIDDIR(N:N) = CHAR(0)
      N = LENC(HELPDIR)+1
      HELPDIR(N:N) = CHAR(0)
      BLANK(1) = X_BVAL
      BLANK(2) = X_EVAL
      CALL MVIEW (X_FILE,MEMORY(IPX),DIM, USER_COORD,
     $UIDDIR,HELPDIR,PIDS,BLANK)
      CALL GDF_FRIS (X_ISLO,ERROR)
*
      END
*<FF>
      SUBROUTINE USER_COORD(IX,JY,KZ,U,CHAINX,CHAINY,CHAINZ,CHAINV)
      CHARACTER*80 CHAINX,CHAINY,CHAINZ,CHAINV
      INTEGER IX,JY,KZ
      REAL *8 U(3)
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'mview.inc'
      REAL*8 UX,VY,WZ,AUX,AVY
      REAL*4 R, VALUE
      INTEGER*4 POINTER
      INTEGER LENC
      CHARACTER*80 CHAIN(3)
*
      DATA CHAIN/' ',' ',' '/
*
* Get value first
      IPX = POINTER(X_ADDR,MEMORY)
      R = VALUE (MEMORY(IPX),DIM(1),DIM(2),DIM(3),
     $IX,JY,KZ)
      WRITE(CHAINV,'(1PG11.3)') R
      CHAINV(13:) = X_UNIT
*
* Convert to original pixels
      IX = IX+OX-1
      JY = JY+OY-1
      KZ = KZ+OZ-1
      U(1) = (DFLOAT(IX)-X_REF1)*X_INC1+X_VAL1
      U(2) = (DFLOAT(JY)-X_REF2)*X_INC2+X_VAL2
      U(3) = (DFLOAT(KZ)-X_REF3)*X_INC3+X_VAL3
*
* Compute the various permutations...
      UX = U(X_XAXI)
      VY = U(X_YAXI)
      WZ = U(X_FAXI)
*
      CALL REL_TO_ABS (UX, VY, AUX, AVY, 1)
*
      CALL SEXAG(CHAIN(X_XAXI),AUX,24)
      CALL SEXAG(CHAIN(X_YAXI),AVY,360)
      WRITE(CHAIN(X_FAXI),'(F12.3)') WZ
*
* Inverse permutation
      CHAINX = CHAIN(1)(1:LENC(CHAIN(1)))//CHAR(0)
      CHAINY = CHAIN(2)(1:LENC(CHAIN(2)))//CHAR(0)
      CHAINZ = CHAIN(3)(1:LENC(CHAIN(3)))//CHAR(0)
      END
*<FF>
      FUNCTION VALUE(Z,NX,NY,NZ,IX,IY,IZ)
      INTEGER NX,NY,NZ,IX,IY,IZ
      REAL VALUE,Z(NX,NY,NZ)
*
      VALUE = Z(IX,IY,IZ)
      END
*<FF>
      SUBROUTINE SPECTRE(NP,POLY,S1)
      INTEGER NP
      REAL POLY(2,NP), S1(*)
      INTEGER S0(4096)
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'mview.inc'
*
      INTEGER*4 POINTER, ADDR, IPA
      INTEGER MP, IER, SIC_GETVM
      INTEGER IMIN,IMAX,JMIN,JMAX,I,J,K
      REAL*8 X,Y,BOUND(5)
      REAL*4 VAL, VALUE, BB, EB
      LOGICAL GR8_IN
*
      DO I=1,DIM(3)
         S0(I) = 0
         S1(I) = 0.
      ENDDO
      MP = NP+1
      IER = SIC_GETVM(8*MP,ADDR)
      IPA = POINTER(ADDR,MEMORY)
      CALL DEFINE_POLY(POLY,NP,MP,MEMORY(IPA),BOUND)
*
* Avoid exploring all the Map by finding IMIN,IMAX,JMIN,JMAX
      IMIN = MAX (1, INT(BOUND(2)) )
      IMAX = MIN (DIM(1),INT(BOUND(3))+1 )
      JMIN = MAX (1, INT(BOUND(4)) )
      JMAX = MIN (DIM(2),INT(BOUND(5))+1 )
*
* Now explore a reasonable part of the map
      IPX = POINTER(X_ADDR,MEMORY)
      IF (X_EVAL.LT.0.0) THEN
         DO J=JMIN,JMAX
            DO I=IMIN,IMAX
               X = DFLOAT(I)
               Y = DFLOAT(J)
               IF (GR8_IN(X,Y,MP,NP,MEMORY(IPA),BOUND)) THEN
                  DO K=1,DIM(3)
                     S0(K) = S0(K) + 1
                     S1(K) = S1(K) + VALUE (MEMORY(IPX),
     $               DIM(1),DIM(2),DIM(3),
     $               I,J,K)
                  ENDDO
               ENDIF
            ENDDO
         ENDDO
      ELSE
         EB = X_EVAL
         BB = X_BVAL
         DO J=JMIN,JMAX
            DO I=IMIN,IMAX
               X = DFLOAT(I)
               Y = DFLOAT(J)
               IF (GR8_IN(X,Y,MP,NP,MEMORY(IPA),BOUND)) THEN
                  DO K=1,DIM(3)
                     VAL = VALUE (MEMORY(IPX),
     $               DIM(1),DIM(2),DIM(3),
     $               I,J,K)
                     IF (ABS(VAL-BB).GT.EB) THEN
                        S0(K) = S0(K) + 1
                        S1(K) = S1(K) + VAL
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
         ENDDO
      ENDIF
*
      X = ABS(X_INC1*X_INC2)
      DO K=1,DIM(3)
         IF (S0(K).NE.0) THEN
            S1(K) = S1(K) / FLOAT(S0(K)) * X
         ENDIF
      ENDDO
      END
*<FF>
      SUBROUTINE DEFINE_POLY (POLY,NP,MP,GONS,BOUND)
      INTEGER NP,MP
      REAL*4 POLY(2,NP)
      REAL*8 GONS(MP,4),BOUND(5)
*
      INTEGER I,NGON
      REAL*8 XGON1,XGON2,YGON1,YGON2
*
* Setup the polygon data structure
      NGON = NP
      DO I=1,NGON
         GONS(I,1) = POLY(1,I)
         GONS(I,2) = POLY(2,I)
      ENDDO
*
      GONS(NGON+1,1) = GONS(1,1)
      GONS(NGON+1,2) = GONS(1,2)
      XGON1 = GONS(1,1)
      XGON2 = GONS(1,1)
      YGON1 = GONS(1,2)
      YGON2 = GONS(1,2)
      DO I=1,NGON
         GONS(I,3) = GONS(I+1,1)-GONS(I,1)
         IF (GONS(I+1,1).LT.XGON1) THEN
            XGON1 = GONS(I+1,1)
         ELSEIF (GONS(I+1,1).GT.XGON2) THEN
            XGON2 = GONS(I+1,1)
         ENDIF
         GONS(I,4) = GONS(I+1,2)-GONS(I,2)
         IF (GONS(I+1,2).LT.YGON1) THEN
            YGON1 = GONS(I+1,2)
         ELSEIF (GONS(I+1,2).GT.YGON2) THEN
            YGON2 = GONS(I+1,2)
         ENDIF
      ENDDO
      BOUND(1) = XGON1-0.01*(XGON2-XGON1)
      BOUND(2) = XGON1
      BOUND(3) = XGON2
      BOUND(4) = YGON1
      BOUND(5) = YGON2
      END
*<FF>
      FUNCTION GR8_IN (X,Y,MGON,NGON,GONS,BOUND)
C----------------------------------------------------------------------
C GREG	Internal routine
C	Find if a point is within a n-gon
C Arguments :
C     X	R*8	Coordinate of point to test	        Input
C     Y	R*8					        Input
C     MGON	I	Dimension of GONS		Input
C     NGON	I	Current number of summits	Input
C     GONS	R*8(*)	Description of polygon		Input
C     BOUND	R*8(5)	Polygon boundary   	        Input
C     GR8_IN	L	Point is within ? 		Output
C No subroutine referenced
C----------------------------------------------------------------------
      INTEGER INDEX,I,MGON,NGON
      REAL*8 X,Y,GONS(MGON,4),BOUND(5)
      LOGICAL GR8_IN
*
      REAL*8 XOUT, D, XX
      LOGICAL YIN
*
* Quick check.
      IF (X.LT.BOUND(2).OR.X.GT.BOUND(3) .OR.
     $Y.LT.BOUND(4).OR.Y.GT.BOUND(5)) THEN
         GR8_IN = .FALSE.
         RETURN
      ENDIF
      GR8_IN = .TRUE.
*
* Count intersection. Note that XOUT is always less than X
      INDEX = 0
      XOUT = BOUND(2)-0.01*(BOUND(3)-BOUND(2))
      DO I=1,NGON
*
         IF (X.EQ.GONS(I,1) .AND. Y.EQ.GONS(I,2)) THEN
            RETURN                       ! Summit is in
         ENDIF
*
* Check Y range
         IF (Y.LE.GONS(I,2) .AND. Y.GE.GONS(I+1,2)) THEN
            YIN = .TRUE.
         ELSEIF (Y.GE.GONS(I,2) .AND. Y.LE.GONS(I+1,2)) THEN
            YIN = .TRUE.
         ELSE
            YIN = .FALSE.
         ENDIF
*
* If Y is OK, check X
         IF (YIN) THEN
*
* Polygon side is normal
            IF (GONS(I,3).NE.0) THEN
               D = GONS(I,4)/GONS(I,3)
* Lines intersect
               IF (D.NE.0D0) THEN
                  XX =  (Y-GONS(I,2))/D + GONS(I,1)
                  IF (XX.LT.X) THEN
                     INDEX = INDEX+1
                  ELSEIF (XX.EQ.X) THEN
                     RETURN
                  ENDIF
* Lines have same slope (i.e. Horizontal)
               ELSE
                  IF (X.GE.GONS(I,1) .AND. X.LE.GONS(I+1,1) ) THEN
                     RETURN              ! On edge
                  ELSEIF (X.LE.GONS(I,1) .AND. X.GE.GONS(I+1,1) ) THEN
                     RETURN              ! On edge
                  ENDIF
               ENDIF
            ELSE
*
* Polygon side is vertical
               IF (GONS(I,1).LT.X) THEN
                  INDEX = INDEX+1        ! Crosses
               ELSEIF (GONS(I,1).EQ.X) THEN
                  RETURN                 ! On edge
               ENDIF
            ENDIF
         ENDIF
      ENDDO
*
      GR8_IN = MOD(INDEX,2).EQ.1
      END
