* Last processed by NICE on 12-Jun-2000 15:50:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE GR8_PROJEC(ICODE,PROJ)
      INCLUDE 'error.inc'
      INCLUDE 'inc:errcod.inc'
      REAL*8 PROJ(3)
      INTEGER ICODE
*
      IF (ICODE.LT.0 .OR. ICODE.GT.7) THEN
         CALL GAGOUT ('F-GR8_PROJEC,  Unvalid projection code ')
         CALL SYSEXI(FATALE)
      ENDIF
      CALL SETPRO(PROJ(1),PROJ(2),PROJ(3),ICODE)
      END
*<FF>
      SUBROUTINE EQ1950_GAL(A,D,LII,BII,N)
C----------------------------------------------------------------------
C ASTRO	Utility routine
C	Converts (A,D) (1950.0) in LII,BII
C	Units are radians
C----------------------------------------------------------------------
      REAL*8 A(1),D(1),LII(1),BII(1)
      INCLUDE 'inc:pi.inc'
      REAL*8 AA,COSA,SINA,COSDEC,SINDEC,EXPR,C626,S626
      PARAMETER (C626=0.4601997847838517D0, S626=0.8878153851364013D0)
      INTEGER I,N
*
* Formulaes are in degrees (LANG K.R. page 504)
*
*	SIN(Bii) 	     = SIN(D)COS(62.6) - COS(D)SIN(A-282.25)SIN(62.6)
*	COS(Bii)*COS(Lii-33) = COS(D)COS(A-282.25)
*	COS(Bii)*SIN(Lii-33) = COS(D)SIN(A-282.25)COS(62.6) + SIN(D)SIN(62.6)
*
      DO I=1,N
         AA = A(I)-282.25D0*PI/180.D0
         COSA = COS(AA)
         SINA = SIN(AA)
         COSDEC = COS(D(I))
         SINDEC = SIN(D(I))
         EXPR = COSDEC*SINA
         BII(I) = ASIN (SINDEC*C626 - EXPR*S626)
         LII(I) = ATAN2( (EXPR*C626+SINDEC*S626) , COSDEC*COSA )
     $   + 33.D0*PI/180.D0
      ENDDO
      END
*<FF>
      SUBROUTINE GAL_EQ1950(LII,BII,A,D,N)
C----------------------------------------------------------------------
C ASTRO	Utility routine
C	Converts Lii,Bii to (A,D) (1950.0)
C	Units are radians
C----------------------------------------------------------------------
      REAL*8 A(1),D(1),LII(1),BII(1)
      INCLUDE 'inc:pi.inc'
      REAL*8 LL,SINB,COSB,SINL,COSL,EXPR,C626,S626
      PARAMETER (C626=0.4601997847838517D0, S626=0.8878153851364013D0)
      INTEGER I,N
*
* The conversion formulae are
*	COS(D)SIN(A-282.25) = COS(Bii)SIN(Lii-33)COS(62.6) - SIN(Bii)SIN(62.6)
*	COS(D)COS(A-282.25) = COS(Bii)COS(Lii-33)
*	SIN(D)		    = COS(Bii)SIN(Lii-33)SIN(62.6) + SIN(Bii)COS(62.6)
*
      DO I = 1,N
         LL   = LII(I) - 33.D0*PI/180.D0
         SINB = SIN(BII(I))
         COSB = COS(BII(I))
         SINL = SIN(LL)
         COSL = COS(LL)
         EXPR = COSB*SINL
         D(I) = ASIN ( EXPR*S626 + SINB*C626 )
         A(I) = ATAN2 ( EXPR*C626 - SINB*S626 , COSB*COSL )
     $   + 282.25D0*PI/180.D0
      ENDDO
      END
*<FF>
      SUBROUTINE REL_TO_ABS(X,Y,A,D,N)
C----------------------------------------------------------------------
C MAP	Internal routine
C	Convert a set of plane position to spherical coordinates
C	according to the current projection formulae
C
C	This routine has been carefully programmed to allow
C	equivalencing (X,A) and (Y,D) arrays. This feature is
C	used in CONVERT.
C
C Arguments :
C	X	R*8 (*)	Array of offsets in X		Input
C	Y	R*8 (*)	Array of offsets in Y		Input
C	A	R*8 (*)	Array of R.A. (radians)		Output
C	D	R*8 (*)	Array of declinations		Output
C	N	I	Size of arrays			Input
C----------------------------------------------------------------------
      INCLUDE 'projec.inc'
      INCLUDE 'inc:pi.inc'
      REAL*8 PRECISION
      PARAMETER (PRECISION=1.0D-10)
      REAL*8 X(*),Y(*),A(*),D(*)
      REAL*8 R,P,SINR,COSR,SINP,COSP,SINDEC,DA,CPANG,SPANG
      INTEGER I,N
      LOGICAL CHANGESIGN
      REAL*8 ARCSIN,ARCCOS
*
      GOTO (001,100,200,300,400,500,600,700) P_TYPE+1
      RETURN
*
* No projection
001   CONTINUE
      IF (ABS(PANGLE).LE.PRECISION) THEN
         DO I = 1,N
            A(I) = A0+X(I)
            D(I) = D0+Y(I)
         ENDDO
      ELSE
         CPANG = COS(PANGLE)
         SPANG = SIN(PANGLE)
         DO I = 1,N
            R = A0+X(I)*CPANG+Y(I)*SPANG
            P = D0-X(I)*SPANG+Y(I)*CPANG
            A(I) = R
            D(I) = P
         ENDDO
      ENDIF
      RETURN
*
* (A,D) = f(A0,D0,R,P)
*	D = ASIN (SIN(D0)*COS(R) + COS(D0)*SIN(R)*COS(P))
*	A = A0 + ATAN2( SIN(R)*SIN(P)*SIN(D0) ,
*	1	 (SIN(D)*COS(D0)-SIN(R)*COS(P)) )
*	ou A = A0 + ASIN ( SIN(R)*SIN(P)/COS(D) )
*
* Gnomonic projection
100   CONTINUE
      IF (SIND0.GT.PRECISION) THEN
         DO I = 1,N
            R = ATAN ( SQRT(X(I)**2+Y(I)**2) )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               IF (Y(I) .LT. NPOLE) THEN
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 + ARCSIN ( SINR*SINP/COS(D(I)) )
               ELSE
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 - ARCSIN ( SINR*SINP/COS(D(I)) ) + PI
               ENDIF
            ELSE
               A(I) = A0
               D(I) = D0
            ENDIF
         ENDDO
      ELSEIF (SIND0.LT.-PRECISION) THEN
         DO I = 1,N
            R = ATAN ( SQRT(X(I)**2+Y(I)**2) )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               IF (Y(I) .GT. SPOLE) THEN
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 + ARCSIN ( SINR*SINP/COS(D(I)) )
               ELSE
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 - ARCSIN ( SINR*SINP/COS(D(I)) ) + PI
               ENDIF
            ELSE
               A(I) = A0
               D(I) = D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            R = ATAN ( SQRT(X(I)**2+Y(I)**2) )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ARCSIN ( SINR*SINP/COS(D(I)) )
            ELSE
               A(I) = A0
               D(I) = D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Orthographic (Dixon)
200   CONTINUE
*
      IF (SIND0.GT.PRECISION) THEN
         DO I = 1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.1.D0) THEN
               SINR = SQRT(R)
               R = ARCSIN ( SINR )
               IF (SINR.GT.PRECISION) THEN
                  P = ATAN2 ( X(I),Y(I) ) - PANGLE
                  COSR = COS(R)
                  SINP = SIN(P)
                  COSP = COS(P)
                  SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $            SINDEC*COSD0-SINR*COSP )
                  IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
               ELSE
                  A(I) = A0
                  D(I) = D0
               ENDIF
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ELSEIF (SIND0.LT.-PRECISION) THEN
         DO I = 1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.1.D0) THEN
               SINR = SQRT(R)
               R = ARCSIN ( SINR )
               IF (SINR.GT.PRECISION) THEN
                  P = ATAN2 ( X(I),Y(I) ) - PANGLE
                  COSR = COS(R)
                  SINP = SIN(P)
                  COSP = COS(P)
                  SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $            SINDEC*COSD0-SINR*COSP ) + PI
                  IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
               ELSE
                  A(I) = A0
                  D(I) = D0
               ENDIF
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.1.D0) THEN
               SINR = SQRT(R)
               R = ARCSIN ( SINR )
               IF (SINR.GT.PRECISION) THEN
                  P = ATAN2 ( X(I),Y(I) ) - PANGLE
                  COSR = COS(R)
                  SINP = SIN(P)
                  COSP = COS(P)
                  SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
                  D(I) = ARCSIN (SINDEC)
                  A(I) = A0 + ATAN2( SINR*SINP,
     $            COSR + SIND0*SINR*COSP)
               ELSE
                  A(I) = A0
                  D(I) = D0
               ENDIF
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Azimuthal  ou Schmidt
300   CONTINUE
      IF (SIND0.GT.PRECISION) THEN
         DO I = 1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.PI-PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP )
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSEIF (R.LT.PI+PRECISION) THEN
               A(I) = A0
               D(I) = -D0
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ELSEIF (SIND0.LT.-PRECISION) THEN
         DO I = 1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.PI-PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP ) + PI
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSEIF (R.LT.PI+PRECISION) THEN
               A(I) = A0
               D(I) = -D0
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I = 1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.PI-PRECISION) THEN
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP,
     $         COSR + SIND0*SINR*COSP)
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSEIF (R.LT.PI+PRECISION) THEN
               A(I) = A0
               D(I) = -D0
            ELSE
               A(I) = 0.D0
               D(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Stereographic
400   CONTINUE
      IF (SIND0.GT.PRECISION) THEN
         DO I = 1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSE
               R = 2.D0 * ATAN (R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP )
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ENDIF
         ENDDO
      ELSEIF (SIND0.LT.-PRECISION) THEN
         DO I = 1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSE
               R = 2.D0 * ATAN (R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP ) + PI
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            R = SQRT(X(I)**2+Y(I)**2)
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSE
               R = 2.D0 * ATAN (R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               IF (0.5D0*PI-ABS(D(I)).GT.PRECISION) THEN
                  A(I) = A0 + ATAN2( SINR*SINP,
     $            COSR + SIND0*SINR*COSP)
               ELSE
                  A(I) = A0
               ENDIF
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Lambert equal area
500   CONTINUE
*     	U = 2*sin(R)/SQRT(2*(1+cos(R))
*	R = Acos (1-0.5*U**2)
      IF (SIND0.GT.PRECISION) THEN
         DO I = 1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.4D0) THEN
               R = ACOS (1D0-0.5D0*R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP )
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSE
               A(I) = 0.0
               D(I) = 0.0
            ENDIF
         ENDDO
      ELSEIF (SIND0.LT.-PRECISION) THEN
         DO I = 1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.4D0) THEN
               R = ACOS (1D0-0.5D0*R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               A(I) = A0 + ATAN2( SINR*SINP*SIND0,
     $         SINDEC*COSD0-SINR*COSP ) + PI
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSE
               A(I) = 0.0
               D(I) = 0.0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            R = X(I)**2+Y(I)**2
            IF (R.LE.PRECISION) THEN
               A(I) = A0
               D(I) = D0
            ELSEIF (R.LE.4D0) THEN
               R = ACOS (1.D0-0.5D0*R)
               P = ATAN2 ( X(I),Y(I) ) - PANGLE
               SINR = SIN(R)
               COSR = COS(R)
               SINP = SIN(P)
               COSP = COS(P)
               SINDEC = (SIND0*COSR + COSD0*SINR*COSP)
               D(I) = ARCSIN (SINDEC)
               IF (0.5D0*PI-ABS(D(I)).GT.PRECISION) THEN
                  A(I) = A0 + ATAN2( SINR*SINP,
     $            COSR + SIND0*SINR*COSP)
               ELSE
                  A(I) = A0
               ENDIF
               IF (A(I).GT.2.D0*PI) A(I) = A(I)-2.D0*PI
            ELSE
               A(I) = 0.0
               D(I) = 0.0
            ENDIF
         ENDDO
      ENDIF
      RETURN
* Aitoff
600   DO I = 1,N
         CHANGESIGN = Y(I).LT.0.D0
         SINR = 0.5D0*SQRT(0.25D0*X(I)**2+Y(I)**2)
         IF (SINR.GT.1.D0) THEN
            A(I) = A0
            D(I) = 0.5D0*PI
         ELSEIF (SINR.GT.PRECISION) THEN
            R = 2.D0 * ARCSIN(SINR)
            P = ATAN2 ( X(I),2.D0*Y(I) )
            COSR = COS(R)
            DA = ATAN2 (SIN(P)*SIN(R),COSR)
            D(I) = ARCCOS ( COSR/COS(DA) )
            A(I) = A0 + 2.D0*DA
         ELSE
            A(I) = A0                    !	+ 2.D0*ARCSIN(SINR)
            D(I) = 0.D0
         ENDIF
         IF (CHANGESIGN) D(I) = -D(I)
      ENDDO
      RETURN
*
* Classical Single-Dish mapping
700   CONTINUE
      DO I=1,N
         IF (Y(I).LE.SPOLE) THEN
            A(I) = A0
            D(I) = -PI*0.5D0
         ELSEIF (Y(I).GE.NPOLE) THEN
            A(I) = A0
            D(I) = PI*0.5D0
         ELSE
            P = D0 + Y(I)
            D(I) = P
            R = X(I)/COS(P)
            IF (R.LE.-PI) THEN
               A(I) = A0-PI
            ELSEIF (R.GE.PI) THEN
               A(I) = A0+PI
            ELSE
               A(I) = A0 + R
            ENDIF
         ENDIF
      ENDDO
      RETURN
      END
*<FF>
      SUBROUTINE ABS_TO_REL(A,D,X,Y,N)
C----------------------------------------------------------------------
C MAP	Internal routine
C	Convert a set of plane position to spherical coordinates
C	according to the current projection formulae
C
C	This routine has been carefully programmed to allow
C	equivalencing (A,X) and (D,Y) arrays. This feature is
C	used in CONVERT.
C
C Arguments :
C	A	R*8 (*)	Array of R.A. (radians)		Input
C	D	R*8 (*)	Array of declinations		Input
C	X	R*8 (*)	Array of offsets in X		Output
C	Y	R*8 (*)	Array of offsets in Y		Output
C	N	I	Size of arrays			Input
C----------------------------------------------------------------------
      INCLUDE 'projec.inc'
      INCLUDE 'inc:pi.inc'
      REAL*8 PRECISION
      PARAMETER (PRECISION=1.0D-10)
      REAL*8 X(*),Y(*),A(*),D(*)
      REAL*8 R,P,COSDEC,SINDEC,COSA,SINA,EXPR,DA
      INTEGER I,N
      LOGICAL CHANGESIGN
      REAL*8 ARCSIN,ARCCOS
*
      GOTO (001,100,200,300,400,500,600,700) P_TYPE+1
      RETURN
*
*
* (R,P) = f(A0,D0,A,D)
*	R = ACOS ( (SIN(D) - COS(D0) *
*	1 ( SIN(D)*COS(D0)-COS(D)*SIN(D0)*COS(A-A0) ) ) / SIN(D0) )
*	P = ATAN2 ( COS(D)*SIN(A-A0) ,
*	1 ( SIN(D)*COS(D0) - COS(D)*SIN(D0)*COS(A-A0) ) )
*
*
*
* Case D0 very small
*	R = ACOS (COS(D)*COS(A-A0))
*	P = A-A0
*
* No projection
001   CONTINUE
      DO I = 1, N
         X(I) = A(I)-A0
         Y(I) = D(I)-D0
      ENDDO
      RETURN
*
* Gnomonic
100   CONTINUE
      IF (ABS(SIND0).GE.PRECISION) THEN
         DO I = 1,N
            SINDEC = SIN(D(I))
            COSDEC = COS(D(I))
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            EXPR = SINDEC*COSD0-COSDEC*SIND0*COSA
            R = ARCCOS ( (SINDEC - COSD0*EXPR)/SIND0 )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA , EXPR ) + PANGLE
               EXPR = TAN(R)
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            COSDEC = COS(D(I))
            SINDEC = SIN(D(I))
            R = ARCCOS ( COSDEC*COSA )
            IF (R.GE.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA,SINDEC ) + PANGLE
               EXPR = TAN(R)
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Dixon ou Orthographic
200   CONTINUE
      IF (ABS(SIND0).GE.PRECISION) THEN
         DO I = 1,N
            SINDEC = SIN(D(I))
            COSDEC = COS(D(I))
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            EXPR = SINDEC*COSD0-COSDEC*SIND0*COSA
            R = ARCCOS ( (SINDEC - COSD0*EXPR)/SIND0 )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA , EXPR ) + PANGLE
               EXPR = SIN(R)
*!!		    if (r.ge.pi/2) r=pi/2
               X(I) = EXPR*SIN(P)
               Y(I) = EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            COSDEC = COS(D(I))
            SINDEC = SIN(D(I))
            R = ARCCOS ( COSDEC*COSA )
            IF (R.GE.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA,SINDEC ) + PANGLE
               EXPR = SIN(R)
*!!		    if (r.ge.pi/2) r=pi/2
               X(I) = EXPR*SIN(P)
               Y(I) = EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Azimuthal or Schmidt
300   CONTINUE
      IF (ABS(SIND0).GE.PRECISION) THEN
         DO I = 1,N
            SINDEC = SIN(D(I))
            COSDEC = COS(D(I))
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            EXPR = SINDEC*COSD0-COSDEC*SIND0*COSA
            R = ARCCOS ( (SINDEC - COSD0*EXPR)/SIND0 )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA , EXPR ) + PANGLE
               X(I) = R*SIN(P)
               Y(I) = R*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            COSDEC = COS(D(I))
            SINDEC = SIN(D(I))
            R = ARCCOS ( COSDEC*COSA )
            IF (R.GE.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA,SINDEC ) + PANGLE
               X(I) = R*SIN(P)
               Y(I) = R*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Stereographic
400   CONTINUE
      IF (ABS(SIND0).GE.PRECISION) THEN
         DO I = 1,N
            SINDEC = SIN(D(I))
            COSDEC = COS(D(I))
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            EXPR = SINDEC*COSD0-COSDEC*SIND0*COSA
            R = ARCCOS ( (SINDEC - COSD0*EXPR)/SIND0 )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA , EXPR ) + PANGLE
               EXPR = TAN(R*0.5D0)
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            COSDEC = COS(D(I))
            SINDEC = SIN(D(I))
            R = ARCCOS ( COSDEC*COSA )
            IF (R.GE.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA,SINDEC ) + PANGLE
               EXPR = TAN(R*0.5D0)
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Lambert equal area
500   CONTINUE
      IF (ABS(SIND0).GE.PRECISION) THEN
         DO I = 1,N
            SINDEC = SIN(D(I))
            COSDEC = COS(D(I))
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            EXPR = SINDEC*COSD0-COSDEC*SIND0*COSA
            R = ARCCOS ( (SINDEC - COSD0*EXPR)/SIND0 )
            IF (R.GT.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA , EXPR ) + PANGLE
               IF (R.GE.PI-PRECISION) THEN
                  EXPR = 2.0
               ELSE
                  EXPR = 2.D0*SIN(R)/SQRT(2.D0*(1.D0+COS(R)))
               ENDIF
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ELSE
         DO I=1,N
            SINA = SIN(A(I)-A0)
            COSA = COS(A(I)-A0)
            COSDEC = COS(D(I))
            SINDEC = SIN(D(I))
            R = ARCCOS ( COSDEC*COSA )
            IF (R.GE.PRECISION) THEN
               P = ATAN2 ( COSDEC*SINA,SINDEC ) + PANGLE
               IF (R.GE.PI-PRECISION) THEN
                  EXPR = 2.0
               ELSE
                  EXPR = 2.D0*SIN(R)/SQRT(2.D0*(1.D0+COS(R)))
               ENDIF
               X(I) =  EXPR*SIN(P)
               Y(I) =  EXPR*COS(P)
            ELSE
               X(I) = 0.D0
               Y(I) = 0.D0
            ENDIF
         ENDDO
      ENDIF
      RETURN
*
* Aitoff
600   CONTINUE
      DO I=1,N
         CHANGESIGN = D(I).LT.0.D0
         DA = A(I) - A0
         DO WHILE (DA.LT.-PI)
            DA = DA + 2.D0*PI
         ENDDO
         DO WHILE (DA.GT.PI)
            DA = DA - 2.D0*PI
         ENDDO
         DA = DA*0.5D0
         IF (DA.EQ.0.D0 .AND. D(I).EQ.0.D0) THEN
            X(I)=0.D0
            Y(I)=0.D0
         ELSE
            COSDEC = COS(D(I))
            R = ARCCOS( COSDEC*COS(DA) )
            EXPR = MIN ( MAX(COSDEC*SIN(DA)/SIN(R),-1.D0), 1.D0)
            R = 2.D0*SIN(R*0.5D0)
            X(I) = 2.D0*R*EXPR
            Y(I) = R*COS(ARCSIN(EXPR))
            IF (CHANGESIGN) Y(I)=-Y(I)
         ENDIF
      ENDDO
      RETURN
*
* Classical Single-Dish mapping
700   DO I=1,N
         P = D(I)
         DA = A(I)-A0
         DO WHILE (DA.LT.-PI)
            DA = DA+2.D0*PI
         ENDDO
         DO WHILE (DA.GT.PI)
            DA = DA-2.D0*PI
         ENDDO
         X(I) = DA*COS(P)
         Y(I) = P-D0
      ENDDO
      RETURN
      END
*<FF>
      FUNCTION ARCSIN (X)
      REAL*8 X,ARCSIN,PI
      PARAMETER (PI=3.14159265358979323846D0)
*
      IF (X.LE.-1.D0) THEN
         ARCSIN = -PI*0.5D0
      ELSEIF (X.GE.1.D0) THEN
         ARCSIN = PI*0.5D0
      ELSE
         ARCSIN = ASIN(X)
      ENDIF
      END
*<FF>
      FUNCTION ARCCOS (X)
      REAL*8 X,ARCCOS,PI
      PARAMETER (PI=3.14159265358979323846D0)
*
      IF (X.LE.-1.D0) THEN
         ARCCOS = PI
      ELSEIF (X.GE.1.D0) THEN
         ARCCOS = 0.0D0
      ELSE
         ARCCOS = ACOS(X)
      ENDIF
      END
*<FF>
      SUBROUTINE SETPRO(RA,DEC,ANG,TYPE)
C----------------------------------------------------------------------
C MAP	Internal routine
C	Define a projection
C Arguments
C	RA	R*8	RA of projection center
C	DEC	R*8	Dec of projection center
C	ANG	R*8	Angle of projection line
C	TYPE	I	Projection type
C----------------------------------------------------------------------
      INCLUDE 'inc:pi.inc'
      INCLUDE 'projec.inc'
      REAL*8 RA,DEC,ANG,X
      INTEGER TYPE
*
* Setup necessary constants
      P_TYPE = TYPE
      A0    = RA
      D0    = DEC
      PANGLE= ANG
      IF (P_TYPE.EQ.P_NONE) RETURN
*
      SIND0 = SIN(D0)
      COSD0 = COS(D0)
      SINA0 = SIN(A0)
      COSA0 = COS(A0)
*
* Compute POLE position and values of remarkable points
*
      IF (P_TYPE.EQ.P_GNOMONIC) THEN
         IF (D0.LT.-1.D-30) THEN
            SPOLE = 1.D0/TAN(D0)
            NPOLE = 1.D38
         ELSEIF (D0.GT.1.D-30) THEN
            NPOLE = 1.D0/TAN(D0)
            SPOLE = 1.D38
         ELSE
            NPOLE = 1.D38
            SPOLE = 1.D38
         ENDIF
      ELSEIF (P_TYPE.EQ.P_ORTHO) THEN
         IF (D0.LT.1.D-30) THEN
            SPOLE = -COSD0
            NPOLE = 1.D38
         ELSEIF(D0.GT.1.D-30)THEN
            NPOLE = COSD0
            SPOLE = 1.D38
         ELSE
            NPOLE = 1
            SPOLE = -1
         ENDIF
      ELSEIF (P_TYPE.EQ.P_AZIMUTHAL) THEN
         IF (D0.LT.0.D0) THEN
            SPOLE = -D0-PI*0.5D0
            NPOLE = PI*0.5+D0
         ELSE
            SPOLE = D0-PI*0.5D0
            NPOLE = PI*0.5D0-D0
         ENDIF
      ELSEIF (P_TYPE.EQ.P_STEREO) THEN
         IF (ABS(D0).GT.1.D-30) THEN
            CALL ABS_TO_REL(0.D0,PI*0.5D0,X,NPOLE,1)
            CALL ABS_TO_REL(0.D0,-PI*0.5D0,X,SPOLE,1)
         ELSE
            NPOLE = 1.D0
            SPOLE = -1.D0
         ENDIF
      ELSEIF (P_TYPE.EQ.P_LAMBERT) THEN
*	    IF (ABS(D0).GT.1.D-30) THEN
         CALL ABS_TO_REL(0.D0,PI*0.5D0,X,NPOLE,1)
         CALL ABS_TO_REL(0.D0,-PI*0.5D0,X,SPOLE,1)
*	    ELSE
*		NPOLE = 2.D0
*		SPOLE = -2.D0
*	    ENDIF
      ELSEIF (P_TYPE.EQ.P_AITOFF) THEN
         PANGLE = 0.D0
         NPOLE = SQRT(2.0D0)
         SPOLE = -NPOLE
         D0 = 0.D0
         COSD0 = 1.D0
         SIND0 = 0.D0
      ELSEIF (P_TYPE.EQ.P_RADIO) THEN
         PANGLE = 0.D0
         NPOLE = PI*0.5D0 - D0
         SPOLE = -D0 - PI*0.5D0
      ENDIF
      END
