PROGRAM FIT_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 POSA_X_BEST, POSA_Y_BEST, POSB_X_BEST, POSB_Y_BEST, & & FLUX_A_BEST, FLUX_B_BEST REAL*8 CX_ORIG, CY_ORIG, CX, CY, RMS, RMS_BEST, IMA, MOD, SCL REAL*8 CX_A, CY_A REAL*4, ALLOCATABLE :: IMA_R4(:,:) REAL*8, ALLOCATABLE :: MODA(:,:), MODB(:,:) 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, & & IND_IMA, NIT, KP, NP, NA, IPX_MIN, IPX_MAX, IPY_MIN, IPY_MAX, IUER ! IF ( IARGC() < 1 ) THEN WRITE ( 6, * ) 'Usage: fit_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' ELSE IF ( IND_IMA == 2 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092783_230.538_GHz.fits' ELSE IF ( IND_IMA == 3 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092787_244.0_GHz.fits' ELSE IF ( IND_IMA == 4 ) THEN FINAM = '/astrogeo/4C+35.28/orig_images/ALMB00092763_246.0_GHz.fits' END IF ! FILOUT = '/tmp/3.3' 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 ! NIT = 8 N1 = 21 N2 = 21 N3 = 21 N4 = 21 N5 = 101 N6 = 31 ! POSA_X_MIN = -3.0*MAS__TO__RAD POSA_X_MAX = 3.0*MAS__TO__RAD ! POSA_Y_MIN = -3.0*MAS__TO__RAD POSA_Y_MAX = 3.0*MAS__TO__RAD ! POSB_X_MIN = -40.0*MAS__TO__RAD ! dx = -30.25 POSB_X_MAX = -20.0*MAS__TO__RAD ! POSB_Y_MIN = -35.0*MAS__TO__RAD ! dy = -21.82 POSB_Y_MAX = -15.0*MAS__TO__RAD ! FLUX_A_MIN = 0.050 FLUX_A_MAX = 0.150 ! FLUX_B_MIN = 0.000 FLUX_B_MAX = 0.030 ! 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) ) ALLOCATE ( MODA(MAP%DIM1,MAP%DIM2), MODB(MAP%DIM1,MAP%DIM2) ) SIG_MAJ = MAP%BEAM_MAJ/DSQRT ( 8.0D0*DLOG(2.0D0) ) SIG_MIN = MAP%BEAM_MIN/DSQRT ( 8.0D0*DLOG(2.0D0) ) DO 410 J1=1,NIT RMS_BEST = 1.D8 DO 420 J2=1,N1 write ( 6, * ) ' J1= ', J1, ' J2= ', J2 POSA_X = POSA_X_MIN + (J2-1)*(POSA_X_MAX - POSA_X_MIN)/(N1-1) DO 430 J3=1,N2 POSA_Y = POSA_Y_MIN + (J3-1)*(POSA_Y_MAX - POSA_Y_MIN)/(N2-1) DO 440 J4=1,N3 POSB_X = POSB_X_MIN + (J4-1)*(POSB_X_MAX - POSB_X_MIN)/(N3-1) DO 450 J5=1,N4 POSB_Y = POSB_Y_MIN + (J5-1)*(POSB_Y_MAX - POSB_Y_MIN)/(N4-1) DO 460 J6=IPX_MIN,IPX_MAX CX_ORIG = (J6 - (MAP%DIM1/2+1))*MAP%STEP_RA DO 470 J7=IPY_MIN,IPY_MAX CY_ORIG = (J7 - (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) MODA(J7,J6) = DEXP( -( (CX+POSA_X)**2/(2*SIG_MAJ**2) + (CY+POSA_Y)**2/(2*SIG_MIN**2) ) ) MODB(J7,J6) = DEXP( -( (CX+POSB_X)**2/(2*SIG_MAJ**2) + (CY+POSB_Y)**2/(2*SIG_MIN**2) ) ) 470 CONTINUE 460 CONTINUE ! DO 480 J8=1,N5 FLUX_A = FLUX_A_MIN + (J8-1)*(FLUX_A_MAX - FLUX_A_MIN)/(N5-1) DO 490 J9=1,N6 FLUX_B = FLUX_B_MIN + (J9-1)*(FLUX_B_MAX - FLUX_B_MIN)/(N6-1) RMS = 0.0D0 DO 4100 J10=IPX_MIN,IPX_MAX DO 4110 J11=IPY_MIN,IPY_MAX IMA = MAP%IMAGE(J11,J10) MOD = FLUX_A*MODA(J11,J10) + FLUX_B*MODB(J11,J10) RMS = RMS + (IMA - MOD)**2 4110 CONTINUE 4100 CONTINUE ! IF ( RMS < RMS_BEST ) THEN RMS_BEST = RMS FLUX_A_BEST = FLUX_A FLUX_B_BEST = FLUX_B POSB_X_BEST = POSB_X POSB_Y_BEST = POSB_Y POSA_X_BEST = POSA_X POSA_Y_BEST = POSA_Y END IF 490 CONTINUE 480 CONTINUE 450 CONTINUE 440 CONTINUE 430 CONTINUE WRITE ( 6, 110 ) 1.D3*DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ), & & FLUX_A_BEST, FLUX_B_BEST, & & POSA_X_BEST*RAD__TO__MAS, POSA_Y_BEST*RAD__TO__MAS, & & POSB_X_BEST*RAD__TO__MAS, POSB_Y_BEST*RAD__TO__MAS CALL FLUSH ( 6 ) 420 CONTINUE ! POSA_X_INT = (POSA_X_MAX - POSA_X_MIN)/2.0 POSA_X_MIN = POSA_X_BEST - POSA_X_INT/2.0 POSA_X_MAX = POSA_X_BEST + POSA_X_INT/2.0 ! POSA_Y_INT = (POSA_Y_MAX - POSA_Y_MIN)/2.0 POSA_Y_MIN = POSA_Y_BEST - POSA_Y_INT/2.0 POSA_Y_MAX = POSA_Y_BEST + POSA_Y_INT/2.0 ! POSB_X_INT = (POSB_X_MAX - POSB_X_MIN)/2.0 POSB_X_MIN = POSB_X_BEST - POSB_X_INT/2.0 POSB_X_MAX = POSB_X_BEST + POSB_X_INT/2.0 ! POSB_Y_INT = (POSB_Y_MAX - POSB_Y_MIN)/2.0 POSB_Y_MIN = POSB_Y_BEST - POSB_Y_INT/2.0 POSB_Y_MAX = POSB_Y_BEST + POSB_Y_INT/2.0 ! FLUX_A_INT = (FLUX_A_MAX - FLUX_A_MIN)/2.0 FLUX_A_MIN = FLUX_A_BEST - FLUX_A_INT/2.0 FLUX_A_MAX = FLUX_A_BEST + FLUX_A_INT/2.0 ! FLUX_B_INT = (FLUX_B_MAX - FLUX_B_MIN)/2.0 FLUX_B_MIN = FLUX_B_BEST - FLUX_B_INT/2.0 FLUX_B_MAX = FLUX_B_BEST + FLUX_B_INT/2.0 WRITE ( 6, 210 ) J1, 1.D3*DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ) 210 FORMAT ( 'Iteration ', I1, ' Rms: ', F7.4, ' mJy' ) 410 CONTINUE ! IMA_R4 = 0.0 DO 4120 J12=IPX_MIN,IPX_MAX CX_ORIG = (J12 - (MAP%DIM1/2+1))*MAP%STEP_RA DO 4130 J13=IPY_MIN,IPY_MAX CY_ORIG = (J13 - (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) MOD = FLUX_A_BEST*DEXP( -( (CX+CX_A)**2/(2*SIG_MAJ**2) + (CY+CY_A)**2/(2*SIG_MIN**2) ) ) + & & FLUX_B_BEST*DEXP( -( (CX+POSB_X)**2/(2*SIG_MAJ**2) + (CY+POSB_Y)**2/(2*SIG_MIN**2) ) ) IMA_R4(J13-IPY_MIN+1,J12-IPY_MIN+1) = MOD IMA_R4(J13-IPY_MIN+1,J12-IPY_MIN+1) = MAP%IMAGE(J13,J12) - MOD MAP%IMAGE(J13,J12) = MAP%IMAGE(J13,J12) - MOD 4130 CONTINUE 4120 CONTINUE ! ! WRITE ( 6, * ) 'RMS_BEST = ', DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ) ! WRITE ( 6, * ) 'MAP%BEAM_POS_ANG = ',MAP%BEAM_POS_ANG WRITE ( 6, '(A)' ) 'Processed image: '//TRIM(FINAM) WRITE ( 6, 110 ) 1.D3*DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ), & & FLUX_A_BEST, FLUX_B_BEST, & & POSA_X_BEST*RAD__TO__MAS, POSA_Y_BEST*RAD__TO__MAS, & & POSB_X_BEST*RAD__TO__MAS, POSB_Y_BEST*RAD__TO__MAS 110 FORMAT ( 'Rms: ', F7.4, ' mJy', ' Flux: ', F7.4, 1X, F7.4, ' mJy Pos_A: ', F7.3, 1X, F7.3, & & ' Pos_B: ', F7.3, 1X, F7.3, ' mas' ) ! IUER = -1 !! MAP%IMAGE(IPX_MIN+10,IPY_MIN+10) = FLUX_A_BEST ! CALL GEN_MAPPLOT ( MAP, 1, 41, 100.0D0, 6.0D0, 5, & ! & 0, 1, 1, FINAM, FILOUT, IUER ) ! END PROGRAM FIT_4C35_ALMA !#!