* Last processed by NICE on 12-Jun-2000 15:51:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM GRID_PROJECT
C----------------------------------------------------------------------
C GDF	Stand alone program
C	Resample an input CLASS grid to a different projection system
C Subroutines
C	(IMAGE),(SIC)
C----------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:errcod.inc'
      INCLUDE 'inc:memory.inc'
      REAL*8 PI
      PARAMETER (PI=3.141592653589793D0)
*
      CHARACTER*80 NAMEX,NAMEY,CENTER_1,CENTER_2
      LOGICAL ERROR
      INTEGER NKEY,LENC,SIC_GETVM,SYS_CODE,MVOC
      PARAMETER (MVOC=8)
      CHARACTER*16 VOCAB(MVOC),KEY,SYSTEM,PROJ
*
      INTEGER N
      REAL*8 CENTER(2),ANGLE
      INTEGER*4 ADDR, IPX, IPY, X, Y, POINTER
      INTEGER*4 BLC(4), TRC(4)
*
      DATA BLC /4*0/, TRC/4*0/
      DATA VOCAB/'NONE','GNOMONIC','ORTHOGRAPHIC','AZIMUTHAL',
     $'STEREOGRAPHIC','LAMBERT','AITOFF','RADIO'/
*
      CALL GILDAS_OPEN
      CALL GILDAS_CHAR('Y_NAME$',NAMEY)
      CALL GILDAS_CHAR('X_NAME$',NAMEX)
      CALL GILDAS_CHAR('PROJECTION$',PROJ)
      CALL GILDAS_CHAR('SYSTEM$',SYSTEM)
      CALL GILDAS_CHAR('CENTER_1$',CENTER_1)
      CALL GILDAS_CHAR('CENTER_2$',CENTER_2)
      CALL GILDAS_DBLE('ANGLE$',ANGLE,1)
      CALL GILDAS_CLOSE
*
      N = LENC(NAMEY)
      CALL SIC_PARSEF(NAMEY(1:N),Y_FILE,' ','.tab')
      CALL GDF_GEIS (Y_ISLO,ERROR)
      CALL GDF_REIS (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,1000) 'F-GRID_PROJECT,  Cannot read input file'
         GOTO 100
      ENDIF
      IF (Y_FORM.NE.FMT_R4) THEN
         WRITE(6,1000) 'E-GRID_PROJECT,  Only real images supported'
         GOTO 100
      ENDIF
      CALL GDF_READY (Y_ISLO,ERROR)
*
      CALL GDF_CHYX
*
* Now define output
      N = LENC(NAMEX)
      CALL SIC_PARSEF(NAMEX(1:N),X_FILE,' ','.tab')
      IF (PROJ.NE.'UNCHANGED') THEN
         CALL SIC_UPPER(PROJ)
         N = LENC(PROJ)
         IF (N.NE.0) THEN
            CALL SIC_AMBIGS('PROJECTION',PROJ,KEY,NKEY,
     $      VOCAB,MVOC,ERROR)
            IF (ERROR) THEN
               WRITE(6,1000)
     $         'F-GRID_PROJECT,  Invalid projection system'
               GOTO 100
            ENDIF
            X_PTYP = NKEY-1
         ENDIF
      ENDIF
*
* Coordinate system
      IF (SYSTEM.NE.'UNCHANGED') THEN
         CALL SIC_UPPER(SYSTEM)
         IF (SYSTEM.EQ.'GALACTIC') THEN
            X_SYST = SYSTEM
            IF (Y_SYST.EQ.'UNKNOWN') THEN
               WRITE(6,1000) 'F-GRID_PROJECT,  Cannot convert from '
     $         ,' UNKNOWN to GALACTIC'
               GOTO 100
            ELSEIF (Y_SYST.EQ.'EQUATORIAL') THEN
               SYS_CODE = 1
            ELSE
               SYS_CODE = 0
            ENDIF
         ELSEIF (SYSTEM.EQ.'EQUATORIAL') THEN
            X_SYST = SYSTEM
            IF (Y_SYST.EQ.'UNKNOWN') THEN
               WRITE(6,1000) 'F-GRID_PROJECT,  Cannot convert from '
     $         ,' UNKNOWN to EQUATORIAL'
               GOTO 100
            ELSEIF (Y_SYST.EQ.'GALACTIC') THEN
               SYS_CODE = -1
            ELSE
               SYS_CODE = 0
            ENDIF
         ELSEIF (SYSTEM.EQ.'UNKNOWN') THEN
            X_SYST = SYSTEM
            IF (Y_SYST.EQ.'GALACTIC') THEN
               WRITE(6,1000) 'F-GRID_PROJECT,  Cannot convert from '
     $         ,' GALACTIC to UNKNOWN'
               GOTO 100
            ELSEIF (SYSTEM.EQ.'EQUATORIAL') THEN
               WRITE(6,1000) 'F-GRID_PROJECT,  Cannot convert from '
     $         ,' EQUATORIAL to UNKNOWN'
               GOTO 100
            ENDIF
            SYS_CODE = 0
         ENDIF
      ELSE
         SYS_CODE = 0
      ENDIF
*
      IF (SYS_CODE.NE.0) THEN
         WRITE(6,1000)
     $   'E-GRID_PROJECT,  System change not yet supported'
         GOTO 100
      ENDIF
*
      IF (X_SYST.EQ.'EQUATORIAL') THEN
         CALL SIC_DECODE(CENTER_1,CENTER(1),24,ERROR)
      ELSE
         CALL SIC_DECODE(CENTER_1,CENTER(1),360,ERROR)
      ENDIF
      IF (ERROR) GOTO 101
      CALL SIC_DECODE(CENTER_2,CENTER(2),360,ERROR)
      IF (ERROR) GOTO 101
      X_A0 = CENTER(1)
      X_D0 = CENTER(2)
      IF (X_PTYP.EQ.7 .AND. ANGLE.NE.0.0) THEN
         WRITE(6,1000)
     $   'F-GRID_PROJECT,  Angle ignored in Radio projection'
         GOTO 100
      ELSEIF (X_PTYP.EQ.6 .AND. ANGLE.NE.0.0) THEN
         WRITE(6,1000)
     $   'F-GRID_PROJECT,  Angle ignored in Aitoff projection'
         GOTO 100
      ENDIF
      X_PANG = ANGLE*PI/180.0D0
*
* Create output image
      ERROR = .FALSE.
      X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4)
      CALL GDF_GEIS (X_ISLO,ERROR)
      CALL GDF_WRITX(X_ISLO,ERROR)
      CALL GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,1000) 'F-GRID_PROJECT,  Cannot create X file'
         WRITE(6,1000) X_FILE
         GOTO 100
      ENDIF
*
* Get some work space for X,Y
      IF (SIC_GETVM(X_DIM(2)*2,ADDR).NE.1) GOTO 100
      X = POINTER(ADDR,MEMORY)
      IF (SIC_GETVM(X_DIM(2)*2,ADDR).NE.1) GOTO 100
      Y = POINTER(ADDR,MEMORY)
*
      CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
      IF (ERROR) GOTO 100
      CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
      IF (ERROR) GOTO 100
      IPX = POINTER(X_ADDR,MEMORY)
      IPY = POINTER(Y_ADDR,MEMORY)
      CALL R4TOR4 (MEMORY(IPY),MEMORY(IPX),X_DIM(1)*X_DIM(2))
      CALL REGRID (MEMORY(IPY),Y_DIM(1),Y_DIM(2),
     $Y_PTYP,Y_A0,MEMORY(IPX),X_PTYP,X_A0,
     $SYS_CODE,MEMORY(X),MEMORY(Y),1,2)
      CALL GDF_FRIS (X_ISLO,ERROR)
      CALL GDF_FRIS (Y_ISLO,ERROR)
*
      STOP 'S-GRID_PROJECT,  Successful completion'
*
100   CALL SYSEXI (FATALE)
101   WRITE(6,1001)
1001  FORMAT('F-GRID_PROJECT,  Bad syntax in center coordinates')
      CALL SYSEXI (FATALE)
*
C+VMS
c1000  FORMAT(1X,A,A)
C-VMS
C+UNIX
1000  FORMAT(A,A)
C-UNIX
      END
*<FF>
      SUBROUTINE REGRID(A,MC,MV,ATYPE,APROJ,
     $B,BTYPE,BPROJ,CODE,XB,YB,IX,IY)
C----------------------------------------------------------------------
C GREG	Stand alone subroutine
C	Regrid an input map in a given projection to
C	an output map in another projection
C Arguments
C	A	R*4(*)	Input GRID of dimensions MC MV
C	MC,MV   I	Size of A
C	ATYPE	I	Type of projection
C	APROJ	R*8(3)	Projection constants
C			PROJ(1) = A0,	PROJ(2)=D0,	PROJ(3)=Angle
C	B	R*4(*)	Output map of dimensions MC MV
C	BTYPE	I
C	BPROJ	R*8(3)
C	CODE	I	Galactic to Equatorial (-1)
C			or Equatorial to Galactic (1)
C	XB,YB	R*8	Work arrays of dimension MV
C       IX,IY   I       Pointer to columns of X and Y
C----------------------------------------------------------------------
      INTEGER MC,MV,ATYPE,BTYPE,CODE,IX,IY
      REAL*4 A(MC,MV), B(MC,MV)
      REAL*8 XB(MC), YB(MC), APROJ(3), BPROJ(3)
      INTEGER I
*
* Setup Input projection
      CALL GR8_PROJEC (ATYPE,APROJ)
*
* Convert Output projection to Absolute coordinates
      DO I=1,MV
         XB(I) = A(IX,I)
         YB(I) = A(IY,I)
      ENDDO
      CALL REL_TO_ABS(XB,YB,XB,YB,MV)
*      do i=1,10
*         type *,i,xb(i),yb(I)
*      ENDDO
*
* Change coordinate system
      IF (CODE.EQ.-1) THEN
*         CALL GAL_EQU (XB,YB,XB,YB,MV)
      ELSEIF (CODE.EQ.1) THEN
*         CALL EQU_GAL (XB,YB,XB,YB,MV)
      ENDIF
*
* Setup Output projection
      CALL GR8_PROJEC (BTYPE,BPROJ)
*
* Convert Absolute coordinates to input projection
      CALL ABS_TO_REL(XB,YB,XB,YB,MV)
*      do i=1,10
*         type *,i,xb(i),yb(I)
*      ENDDO
      DO I=1,MV
         B(IX,I) = XB(I)
         B(IY,I) = YB(I)
      ENDDO
      END
