PROGRAM PLOT_2COMP_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 FLUX_A, FLUX_B, POSA_X, POSA_Y, POSB_X, POSB_Y, & & SIG_MIN, SIG_MAJ REAL*8 CX_ORIG, CY_ORIG, CX, CY, RMS, IMA, MOD, SCL REAL*8 CX_A, CY_A REAL*4, ALLOCATABLE :: IMA_R4(:,:), MOD_R4(:,:), RES_R4(:,:) CHARACTER STR*128 REAL*8 DX(4), DY(4), DL(4), DX_MEAN, DY_MEAN, DX_RMS, DY_RMS, & & DL_MEAN, DL_RMS, SUM_PXL_ORIG, SUM_PXL_RESID, PXL_MIN, PXL_MAX 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_2comp_4c35_alma ima' CALL EXIT ( 1 ) ELSE CALL GETARG ( 1, STR ) CALL CHIN ( STR, IND_IMA ) END IF DEV = 1 FILOUT = '/tmp/foo.ps' !! IF ( IND_IMA == 1 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092793_228.7_GHz.fits' POSA_X = -0.384*MAS__TO__RAD POSA_Y = 1.205*MAS__TO__RAD POSB_X = -23.750*MAS__TO__RAD POSB_Y = -25.109*MAS__TO__RAD FLUX_A = 0.1026 FLUX_B = 0.0118 RMS = 0.0005802 ELSE IF ( IND_IMA == 2 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092783_230.538_GHz.fits' POSA_X = -0.434*MAS__TO__RAD POSA_Y = 1.228*MAS__TO__RAD POSB_X = -23.008*MAS__TO__RAD POSB_Y = -26.461*MAS__TO__RAD FLUX_A = 0.1056 FLUX_B = 0.0119 RMS = 0.0006228 ELSE IF ( IND_IMA == 3 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092787_244.0_GHz.fits' POSA_X = -0.384*MAS__TO__RAD POSA_Y = 1.184*MAS__TO__RAD POSB_X = -23.906*MAS__TO__RAD POSB_Y = -25.117*MAS__TO__RAD FLUX_A = 0.0991 FLUX_B = 0.0114 RMS = 0.0005554 ELSE IF ( IND_IMA == 4 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092763_246.0_GHz.fits' POSA_X = -0.384*MAS__TO__RAD POSA_Y = 1.242*MAS__TO__RAD POSB_X = -23.477*MAS__TO__RAD POSB_Y = -24.877*MAS__TO__RAD FLUX_A = 0.0984 FLUX_B = 0.0116 RMS = 0.0005576 ! posb_x = (-30.25 + 0.384)*mas__to__rad ! posb_y = (-21.82 - 1.242)*mas__to__rad ! Obs_date: 20231111 Exp_name: ug004o Band: X Total_flux: 0.764 Tot_flux_A: 0.730 Flux_A_peak: 0.655 Flux_B: 0.035 -+ 0.001 Jy Dx: -30.19 -+ 0.01 mas Dy: -21.70 -+ 0.02 mas Dist: 37.18 -+ 0.01 mas Beam: 1.0 x 1.9 bpa: 17.4 A_unr_min: 0.9 A_unr_maj: 0.5 A_bpa: 3.6 & Tb_A: 4.7E+10 Tb_B: >2.4E+09 ! Obs_date: 20151126 Exp_name: 00597.S Band: 6 Total_flux: 0.110 Tot_flux_A: 0.098 Flux_A_peak: 0.098 Flux_B: 0.012 -+ 0.001 Jy Dx: -23.09 -+ 0.50 mas Dy: -26.12 -+ 0.50 mas Dist: 34.86 -+ 0.70 mas Beam: 28.4 x 47.8 bpa: 13.1 A_unr_min: 28.4 A_unr_maj: n/a A_bpa: n/a & Tb_A: >3.E+3 Tb_B: >3.E2 ! 23.09 END IF ! DX(1) = -23.750 + 0.384 DY(1) = -25.109 - 1.205 DX(2) = -23.008 + 0.434 DY(2) = -26.461 - 1.228 DX(3) = -23.906 + 0.384 DY(3) = -25.117 - 1.184 DX(4) = -23.477 + 0.384 DY(4) = -24.877 - 1.242 ! DX_MEAN = (DX(1) + DX(2) + DX(3) + DX(4) )/4 DY_MEAN = (DY(1) + DY(2) + DY(3) + DY(4) )/4 DL_MEAN = ( DSQRT(DX(1)**2 + DY(1)**2) + DSQRT(DX(2)**2 + DY(2)**2) + & & DSQRT(DX(3)**2 + DY(3)**2) + DSQRT(DX(4)**2 + DY(4)**2) )/4 DX_RMS = DSQRT ( ((DX(1) - DX_MEAN)**2 + (DX(2) - DX_MEAN)**2 + (DX(3) - DX_MEAN)**2 + (DX(4) - DX_MEAN)**2)/3 ) DY_RMS = DSQRT ( ((DY(1) - DY_MEAN)**2 + (DY(2) - DY_MEAN)**2 + (DY(3) - DY_MEAN)**2 + (DY(4) - DY_MEAN)**2)/3 ) DL_RMS = DSQRT ( ( (DSQRT(DX(1)**2 + DY(1)**2) - DL_MEAN)**2 + (DSQRT(DX(2)**2 + DY(2)**2) - DL_MEAN)**2 + & & (DSQRT(DX(3)**2 + DY(3)**2) - DL_MEAN)**2 + (DSQRT(DX(4)**2 + DY(4)**2) - DL_MEAN)**2 )/3 ) ! WRITE ( 6, * ) 'DX_MEAN= ', DX_MEAN, ' DX_RMS = ', DX_RMS WRITE ( 6, * ) 'DY_MEAN= ', DY_MEAN, ' DY_RMS = ', DY_RMS WRITE ( 6, * ) 'DL_MEAN= ', DL_MEAN, ' DL_RMS = ', DL_RMS WRITE ( 6, * ) 'DD_RMS = ', DSQRT ( (DX_RMS**2*3 + DY_RMS**2*3)/7 ) ! 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(MAP%BEAM_POS_ANG) + CY_ORIG*DSIN(MAP%BEAM_POS_ANG) CY = -CX_ORIG*DSIN(MAP%BEAM_POS_ANG) + CY_ORIG*DCOS(MAP%BEAM_POS_ANG) IMA_R4(J2,J1) = MAP%IMAGE(J2,J1) MOD_R4(J2,J1) = & & FLUX_A*DEXP( -( (CX+POSA_X)**2/(2*SIG_MAJ**2) + (CY+POSA_Y)**2/(2*SIG_MIN**2) ) ) + & & FLUX_B*DEXP( -( (CX+POSB_X)**2/(2*SIG_MAJ**2) + (CY+POSB_Y)**2/(2*SIG_MIN**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) !! MAP%IMAGE(J4,J3) = IMA_R4(J4,J3) 440 CONTINUE 430 CONTINUE ! !! MAP%IMAGE(IPX_MIN+10,IPY_MIN+10) = FLUX_A_BEST CALL GEN_MAPPLOT ( MAP, DEV, 41, 100.0D0, 5.0D0, 4, & & 0, 1, FINAM, FILOUT, IUER ) ! END PROGRAM PLOT_2COMP_4C35_ALMA !#!