PROGRAM FIT_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_MIN, FLUX_MAX, FLUX_INT REAL*8 POS_X_MIN, POS_X_MAX, POS_X_INT REAL*8 POS_Y_MIN, POS_Y_MAX, POS_Y_INT REAL*8 SIG_X_MIN, SIG_X_MAX, SIG_X_INT REAL*8 SIG_Y_MIN, SIG_Y_MAX, SIG_Y_INT REAL*8 BPA_MIN, BPA_MAX, BPA_INT REAL*8 FLUX, POS_X, POS_Y, SIG_X, SIG_Y, BPA REAL*8 POS_X_BEST, POS_Y_BEST, SIG_X_BEST, SIG_Y_BEST, BPA_BEST, FLUX_BEST REAL*8 CX_ORIG, CY_ORIG, CX, CY, RMS, RMS_BEST, IMA, SCL REAL*8 CX_A, CY_A REAL*4, ALLOCATABLE :: IMA_R4(:,:) REAL*8, ALLOCATABLE :: MOD(:,:) 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 = 6 N1 = 11 N2 = 11 N3 = 21 N4 = 21 N5 = 51 N6 = 51 ! POS_X_MIN = -3.0*MAS__TO__RAD POS_X_MAX = 3.0*MAS__TO__RAD ! POS_Y_MIN = -3.0*MAS__TO__RAD POS_Y_MAX = 3.0*MAS__TO__RAD ! SIG_X_MIN = 10.0D0*MAS__TO__RAD SIG_X_MAX = 40.0D0*MAS__TO__RAD ! SIG_Y_MIN = 10.0D0*MAS__TO__RAD SIG_Y_MAX = 40.0D0*MAS__TO__RAD ! BPA_MIN = 0.0 BPA_MAX = 2.0 ! FLUX_MIN = 0.050 FLUX_MAX = 0.150 ! 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 ( MOD(MAP%DIM1,MAP%DIM2) ) DO 410 J1=1,NIT RMS_BEST = 1.D8 DO 420 J2=1,N1 write ( 6, * ) ' J1= ', J1, ' J2= ', J2 POS_X = POS_X_MIN + (J2-1)*(POS_X_MAX - POS_X_MIN)/(N1-1) DO 430 J3=1,N2 POS_Y = POS_Y_MIN + (J3-1)*(POS_Y_MAX - POS_Y_MIN)/(N2-1) DO 440 J4=1,N3 SIG_X = SIG_X_MIN + (J4-1)*(SIG_X_MAX - SIG_X_MIN)/(N3-1) DO 450 J5=1,N4 SIG_Y = SIG_Y_MIN + (J5-1)*(SIG_Y_MAX - SIG_Y_MIN)/(N4-1) DO 460 J6=1,N5 BPA = BPA_MIN + (J6-1)*(BPA_MAX - BPA_MIN)/(N5-1) DO 470 J7=IPX_MIN,IPX_MAX CX_ORIG = (J7 - (MAP%DIM1/2+1))*MAP%STEP_RA DO 480 J8=IPY_MIN,IPY_MAX CY_ORIG = (J8 - (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) MOD(J8,J7) = DEXP( -( (CX+POS_X)**2/(2*SIG_X**2) + (CY+POS_Y)**2/(2*SIG_Y**2) ) ) 480 CONTINUE 470 CONTINUE ! DO 490 J9=1,N6 FLUX = FLUX_MIN + (J9-1)*(FLUX_MAX - FLUX_MIN)/(N6-1) !! write ( 6, * ) 'f6= ', flux, bpa, sig_x, sig_y, pos_x, pos_y ! %%%%%%%%%%%5 RMS = 0.0D0 DO 4100 J10=IPX_MIN,IPX_MAX DO 4110 J11=IPY_MIN,IPY_MAX IMA = MAP%IMAGE(J11,J10) RMS = RMS + (IMA - FLUX*MOD(J11,J10))**2 4110 CONTINUE 4100 CONTINUE ! IF ( RMS < RMS_BEST ) THEN RMS_BEST = RMS BPA_BEST = BPA FLUX_BEST = FLUX POS_X_BEST = POS_X POS_Y_BEST = POS_Y SIG_X_BEST = SIG_X SIG_Y_BEST = SIG_Y END IF 490 CONTINUE 460 CONTINUE 450 CONTINUE 440 CONTINUE 430 CONTINUE WRITE ( 6, 110 ) 1.D3*DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ), & & FLUX_BEST, & & POS_X_BEST*RAD__TO__MAS, POS_Y_BEST*RAD__TO__MAS, & & SIG_X_BEST*RAD__TO__MAS, SIG_Y_BEST*RAD__TO__MAS, & & BPA_BEST CALL FLUSH ( 6 ) 420 CONTINUE ! POS_X_INT = (POS_X_MAX - POS_X_MIN)/2.0 POS_X_MIN = POS_X_BEST - POS_X_INT/2.0 POS_X_MAX = POS_X_BEST + POS_X_INT/2.0 ! POS_Y_INT = (POS_Y_MAX - POS_Y_MIN)/2.0 POS_Y_MIN = POS_Y_BEST - POS_Y_INT/2.0 POS_Y_MAX = POS_Y_BEST + POS_Y_INT/2.0 ! SIG_X_INT = (SIG_X_MAX - SIG_X_MIN)/2.0 SIG_X_MIN = SIG_X_BEST - SIG_X_INT/2.0 SIG_X_MAX = SIG_X_BEST + SIG_X_INT/2.0 ! SIG_Y_INT = (SIG_Y_MAX - SIG_Y_MIN)/2.0 SIG_Y_MIN = SIG_Y_BEST - SIG_Y_INT/2.0 SIG_Y_MAX = SIG_Y_BEST + SIG_Y_INT/2.0 ! FLUX_INT = (FLUX_MAX - FLUX_MIN)/2.0 FLUX_MIN = FLUX_BEST - FLUX_INT/2.0 FLUX_MAX = FLUX_BEST + FLUX_INT/2.0 ! BPA_INT = (BPA_MAX - BPA_MIN)/2.0 BPA_MIN = BPA_BEST - BPA_INT/2.0 BPA_MAX = BPA_BEST + BPA_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(BPA_BEST) + CY_ORIG*DSIN(BPA_BEST) CY = -CX_ORIG*DSIN(BPA_BEST) + CY_ORIG*DCOS(BPA_BEST) MOD(J13,J2) = FLUX_BEST*DEXP( -( (CX+POS_X_BEST)**2/(2*SIG_X_BEST**2) + & & (CY+POS_Y_BEST)**2/(2*SIG_Y_BEST**2) ) ) IMA_R4(J13-IPY_MIN+1,J12-IPY_MIN+1) = MOD(J13,J12) IMA_R4(J13-IPY_MIN+1,J12-IPY_MIN+1) = MAP%IMAGE(J13,J12) - MOD(J13,J12) MAP%IMAGE(J13,J12) = MAP%IMAGE(J13,J12) - MOD(J13,J12) 4130 CONTINUE 4120 CONTINUE ! WRITE ( 6, '(A)' ) 'Processed image: '//TRIM(FINAM) WRITE ( 6, 110 ) 1.D3*DSQRT( RMS_BEST/(IPX_MAX - IPX_MIN + 1)**2 ), & & FLUX_BEST, & & POS_X_BEST*RAD__TO__MAS, POS_Y_BEST*RAD__TO__MAS, & & SIG_X_BEST*RAD__TO__MAS, SIG_Y_BEST*RAD__TO__MAS, & & BPA_BEST 110 FORMAT ( 'Rms: ', F7.4, ' mJy', ' Flux: ', F7.4, ' mJy ', & & ' Pos: ', F7.3, 1X, F7.3, & & ' Sig: ', F7.3, 1X, F7.3, ' mas', & & ' Bpa: ', F6.3, ' rad' ) ! IUER = -1 CALL GEN_MAPPLOT ( MAP, 1, 41, 100.0D0, 6.0D0, 5, & & 0, 1, 1, FINAM, FILOUT, IUER ) ! END PROGRAM FIT_ELL_4C35_ALMA !#!