PROGRAM    PLOT_ELL_4C35_ALMA
      IMPLICIT   NONE 
      INCLUDE    'astro_constants.i'
      INCLUDE    'sou_map.i'
      LOGICAL*4  FL_IMA, FL_CC
      TYPE     ( SOUMAP__TYPE ) :: MAP
      CHARACTER  FINAM*128, FILOUT*128
      INTEGER*4  MP
      PARAMETER  ( MP = 8192 )
      REAL*8     FLUX_A_MIN, FLUX_A_MAX, FLUX_A_INT
      REAL*8     FLUX_B_MIN, FLUX_B_MAX, FLUX_B_INT
      REAL*8     POSA_X_MIN, POSA_X_MAX, POSA_X_INT
      REAL*8     POSA_Y_MIN, POSA_Y_MAX, POSA_Y_INT
      REAL*8     POSB_X_MIN, POSB_X_MAX, POSB_X_INT
      REAL*8     POSB_Y_MIN, POSB_Y_MAX, POSB_Y_INT
      REAL*8     POS_X, POS_Y, SIG_X, SIG_Y, BPA, FLUX
      REAL*8     CX_ORIG, CY_ORIG, CX, CY, RMS, IMA, MOD, SCL
      REAL*8     CX_A, CY_A, SUM_PXL_ORIG, SUM_PXL_RESID, PXL_MIN, PXL_MAX, &
     &           SIG_MIN, SIG_MAJ
      REAL*4,    ALLOCATABLE :: IMA_R4(:,:), MOD_R4(:,:), RES_R4(:,:)
      CHARACTER  STR*128
      INTEGER*4  N1, N2, N3, N4, N5, N6
      INTEGER*4  J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11, J12, J13, DEV, &
     &           IND_IMA, NIT, KP, NP, NA, IPX_MIN, IPX_MAX, IPY_MIN, IPY_MAX, IUER
!
      IF ( IARGC() < 1 ) THEN
           WRITE ( 6, * ) 'Usage: plot_ell_4c35_alma ima'
           CALL EXIT ( 1 )
         ELSE
           CALL GETARG ( 1, STR )
           CALL CHIN   ( STR, IND_IMA )
      END IF

      IF ( IND_IMA == 1 ) THEN
           FINAM  = '/astrogeo/4C+35.28/orig_images/ALMB00092793_228.7_GHz.fits'
           POS_X = -1.069*MAS__TO__RAD
           POS_Y =  0.488*MAS__TO__RAD
           SIG_X = 20.969*MAS__TO__RAD
           SIG_Y = 12.909*MAS__TO__RAD
           BPA   = 0.319
           FLUX  = 0.1017
           RMS   = 0.0011272
         ELSE IF ( IND_IMA == 2 ) THEN
           FINAM  = '/astrogeo/4C+35.28/orig_images/ALMB00092783_230.538_GHz.fits'
           POS_X =  1.069*MAS__TO__RAD
           POS_Y =  0.506*MAS__TO__RAD
           SIG_X = 13.234*MAS__TO__RAD
           SIG_Y = 21.818*MAS__TO__RAD
           BPA   = 1.840
           FLUX  = 0.1048
           RMS   = 0.0011987
         ELSE IF ( IND_IMA == 3 ) THEN
           FINAM  = '/astrogeo/4C+35.28/orig_images/ALMB00092787_244.0_GHz.fits'
           POS_X = -1.069*MAS__TO__RAD
           POS_Y =  0.488*MAS__TO__RAD
           SIG_X = 20.688*MAS__TO__RAD
           SIG_Y = 12.812*MAS__TO__RAD
           BPA   = 0.328
           FLUX  = 0.0980
           RMS   = 0.0010858
         ELSE IF ( IND_IMA == 4 ) THEN
           FINAM  = '/astrogeo/4C+35.28/orig_images/ALMB00092763_246.0_GHz.fits'
           POS_X = -1.069*MAS__TO__RAD
           POS_Y =  0.488*MAS__TO__RAD
           SIG_X = 20.547*MAS__TO__RAD
           SIG_Y = 12.766*MAS__TO__RAD
           BPA   = 0.329
           FLUX  = 0.0971
           RMS   = 0.0010969
      END IF
!
      FILOUT = '/tmp/foo.ps'
      DEV = 2
!
      FL_CC  = .FALSE.
      FL_IMA = .TRUE.
!
! --- Get the image
!
      IUER = -1
      CALL GET_FITS_MAP ( FINAM, FL_CC, FL_IMA, MAP, IUER )
      IF ( IUER .NE. 0 ) THEN
           CALL ERR_LOG ( 1407, -2, 'FITS_TO_MAP', 'Failure in an '// &
     &         'attempt to read the fits file with image '//FINAM )
           CALL EXIT ( 1 )
      END IF
!
      IPX_MIN = 151 - 20
      IPX_MAX = 151 + 20
      IPY_MIN = 151 - 20
      IPY_MAX = 151 + 20
!
      ALLOCATE ( IMA_R4(1:IPX_MAX-IPX_MIN+1,1:IPY_MAX-IPY_MIN+1), &
     &           MOD_R4(1:IPX_MAX-IPX_MIN+1,1:IPY_MAX-IPY_MIN+1), &
     &           RES_R4(1:IPX_MAX-IPX_MIN+1,1:IPY_MAX-IPY_MIN+1)  )
      IMA_R4 = 0
      MOD_R4 = 0
      RES_R4 = 0
!
      SIG_MAJ = MAP%BEAM_MAJ/DSQRT ( 8.0D0*DLOG(2.0D0) )
      SIG_MIN = MAP%BEAM_MIN/DSQRT ( 8.0D0*DLOG(2.0D0) )
      SUM_PXL_ORIG  = 0.0D0
      SUM_PXL_RESID = 0.0D0
      PXL_MIN =  1.D8
      PXL_MAX = -1.D8
      DO 410 J1=IPX_MIN,IPX_MAX
         CX_ORIG = (J1 - (MAP%DIM1/2+1))*MAP%STEP_RA 
         DO 420 J2=IPY_MIN,IPY_MAX
            CY_ORIG = (J2 - (MAP%DIM2/2+1))*MAP%STEP_DL
            CX =  CX_ORIG*DCOS(BPA) + CY_ORIG*DSIN(BPA) 
            CY = -CX_ORIG*DSIN(BPA) + CY_ORIG*DCOS(BPA) 
            IMA_R4(J2,J1) = MAP%IMAGE(J2,J1) 
            MOD_R4(J2,J1) = FLUX*DEXP( -( (CX+POS_X)**2/(2*SIG_X**2) + (CY+POS_Y)**2/(2*SIG_Y**2) ) )
            RES_R4(J2,J1) = IMA_R4(J2,J1) - MOD_R4(J2,J1) 
            SUM_PXL_ORIG  = SUM_PXL_ORIG  + IMA_R4(J2,J1) 
            SUM_PXL_RESID = SUM_PXL_RESID + RES_R4(J2,J1) 
            IF ( RES_R4(J2,J1) < PXL_MIN ) PXL_MIN = RES_R4(J2,J1) 
            IF ( RES_R4(J2,J1) > PXL_MAX ) PXL_MAX = RES_R4(J2,J1) 
 420     CONTINUE 
 410  CONTINUE 
!
      SCL = DSQRT(SIG_MIN*SIG_MAJ)/MAP%STEP_DL*2*PI__NUM**2
      WRITE  ( 6, 210 ) SUM_PXL_ORIG/SCL, SUM_PXL_RESID/SCL, SUM_PXL_RESID/SUM_PXL_ORIG, &
     &                  PXL_MIN, PXL_MAX
 210  FORMAT ( 'Sum_pxl= ', F7.5, ' Sum_resid: ', F7.5, ' rat= ', F6.4, &
     &         ' PXL_MIN= ', F7.5, ' PXL_MAX= ', F7.5 )
!
      DO 430 J3=IPX_MIN,IPX_MAX
         DO 440 J4=IPY_MIN,IPY_MAX
            MAP%IMAGE(J4,J3) = RES_R4(J4,J3)
 440     CONTINUE 
 430  CONTINUE 
!
      CALL GEN_MAPPLOT ( MAP, DEV, 41, 100.0D0, 5.0D0, 4, &
     &                   0, 1, FINAM, FILOUT, IUER  )
!
      END PROGRAM  PLOT_ELL_4C35_ALMA  !#!