PROGRAM CORE_SHIFT IMPLICIT NONE INCLUDE 'astro_constants.i' INTEGER*4 MP PARAMETER ( MP = 128 ) CHARACTER FILIN*128, BUF(MP)*512, COMP*1, BAND*1, DATE_STR*19, BAND_ARR(MP)*1, & & BANDS*7 REAL*8 DIST(MP), DIST_ERR(MP), DIST_BAND(MP), LOG_FREQ(MP), LOG_DIST(MP), DIST_BAND_ERR(MP), & & FREQ(MP), FREQ_BAND(MP), ERR_MAX, WEI, WEI_SUM, FLOOR_DIST_ERR, & & FREQ_MOD(MP), DIST_MOD(MP), DIST_MOD_ERR(MP), & & DIFF_SX(MP), DIFF_SX_ERR(MP), TIME_SX(MP), & & DIFF_CX(MP), DIFF_CX_ERR(MP), TIME_CX(MP), SEC, & & TIME_K(MP), DIST_K(MP), DIST_K_ERR(MP), & & TIME_S(MP), DIST_S(MP), DIST_S_ERR(MP), & & TIME_C(MP), DIST_C(MP), DIST_C_ERR(MP), & & TIME_X(MP), DIST_X(MP), DIST_X_ERR(MP), DR, SH, SH_ERR, & & DIST_K_W(MP), SH_K, DR_K, DR_K_ERR, MEAN_T_K, & & DIST_X_W(MP), SH_X, DR_X, DR_X_ERR, MEAN_T_X, & & DIST_C_W(MP), SH_C, DR_C, DR_C_ERR, MEAN_T_C, & & DIST_S_W(MP), SH_S, DR_S, DR_S_ERR, MEAN_T_S, & & DX_X(MP), DX_X_ERR(MP), DX_X_W(MP), SH_X_X, DR_X_X, DR_X_X_ERR, MEAN_X_T, & & DY_X(MP), DY_X_ERR(MP), DY_X_W(MP), SH_Y_X, DR_Y_X, DR_Y_X_ERR, MEAN_Y_T, & & TREND_X(MP), TREND_X_ERR(MP), & & DX, DY, DX_ERR, DY_ERR, TRA(MP), TRA_X(MP), & & SH_X_X_ERR, SH_Y_X_ERR, & & ANG, DIR_ANG, DIR_ANG_ERR REAL*8 DIST_SX_BAND(MP,2), DIST_SX_BAND_ERR(MP,2) INTEGER*4 J1, J2, J3, J4, J5, KP, NP, NB, K_SX, K_CX, KS, KC, KX, KK, & & IV(MP), MJD_0, MJD, IUER REAL*8, EXTERNAL :: ATAN_CS ! FILIN = '/astrogeo/4C+35.28/4C+35.28_properties.txt' ERR_MAX = 1.0D0 FLOOR_DIST_ERR = 0.02D0 BANDS = 'LSCXUKQ' ANG = -125.5D0*DEG__TO__RAD ! DATE_STR = '1990.01.01_00:00:00' CALL DATE_TO_TIME ( DATE_STR, MJD_0, SEC, IUER ) ! ! 37.28 -+ 0.05 0.857 ==> 7780ps/arcsec ==> 290.0 -+ 0.4 pc ! IUER = -1 CALL RD_TEXT ( FILIN, MP, BUF, NP, IUER ) IF ( IUER .NE. 0 ) CALL EXIT ( 1 ) ! KP = 0 KS = 0 KC = 0 KX = 0 KK = 0 K_SX = 0 K_CX = 0 DO 410 J1=1,NP IF ( INDEX ( BUF(J1), 'bd152i' ) > 0 ) THEN GOTO 410 END IF KP = KP + 1 IF ( BUF(J1)(46:46) == 'L' ) THEN FREQ(KP) = 1.42D9 ELSE IF ( BUF(J1)(46:46) == 'S' ) THEN FREQ(KP) = 2.30D9 ELSE IF ( BUF(J1)(46:46) == 'C' ) THEN FREQ(KP) = 4.4D9 ELSE IF ( BUF(J1)(46:46) == 'X' ) THEN FREQ(KP) = 8.4D9 ELSE IF ( BUF(J1)(46:46) == 'U' ) THEN FREQ(KP) = 15.2D9 ELSE IF ( BUF(J1)(46:46) == 'K' ) THEN FREQ(KP) = 22.2D9 ELSE IF ( BUF(J1)(46:46) == 'Q' ) THEN FREQ(KP) = 43.2D9 END IF FREQ(KP) = FREQ(KP) + J1*1.D5 READ ( UNIT=BUF(J1)(196:201), FMT='(F6.2)' ) DIST(KP) READ ( UNIT=BUF(J1)(142:147), FMT='(F6.2)' ) DX READ ( UNIT=BUF(J1)(141:156), FMT='(F6.2)' ) DX_ERR READ ( UNIT=BUF(J1)(169:174), FMT='(F6.2)' ) DY READ ( UNIT=BUF(J1)(178:183), FMT='(F6.2)' ) DY_ERR READ ( UNIT=BUF(J1)(206:210), FMT='(F5.2)' ) DIST_ERR(KP) BAND_ARR(KP) = BUF(J1)(46:46) !!! TRA(KP) = DY*DCOS(ANG+P2I) + DX*DSIN(ANG+P2I) !! DIST(KP) = DY*DCOS(ANG -16*deg__to__rad) + DX*DSIN(ANG -16*deg__to__rad) !! DIST(KP) = DY*DCOS(ANG) + DX*DSIN(ANG) IF ( DIST_ERR(KP) > ERR_MAX ) THEN KP = KP - 1 GOTO 410 END IF DATE_STR = BUF(J1)(11:14)//'.'//BUF(J1)(15:16)//'.'//BUF(J1)(17:18)//'_00:00:00' CALL DATE_TO_TIME ( DATE_STR, MJD, SEC, IUER ) IF ( BAND_ARR(KP) == 'S' ) THEN KS = KS + 1 TIME_S(KS) = 1990 + (MJD - MJD_0)/365.25D0 DIST_S(KS) = DIST(KP) DIST_S_ERR(KS) = DIST_ERR(KP) DIST_S_W(KS) = 1.D0/DIST_ERR(KP) IF ( DIST_S_ERR(KS) > 1.8 ) KS = KS - 1 END IF IF ( BAND_ARR(KP) == 'C' ) THEN KC = KC + 1 TIME_C(KC) = 1990 + (MJD - MJD_0)/365.25D0 DIST_C(KC) = DIST(KP) DIST_C_ERR(KC) = DIST_ERR(KP) DIST_C_W(KC) = 1.D0/DIST_ERR(KP) IF ( DIST_C_ERR(KC) > 0.4 ) KC = KC - 1 END IF IF ( BAND_ARR(KP) == 'X' ) THEN KX = KX + 1 TIME_X(KX) = 1990 + (MJD - MJD_0)/365.25D0 DIST_X(KX) = DIST(KP) DIST_X_ERR(KX) = DIST_ERR(KP) DIST_X_W(KX) = 1.D0/DIST_ERR(KP) TRA_X(KX) = DY*DCOS(ANG+P2I) + DX*DSIN(ANG+P2I) DX_X(KX) = DX DX_X_ERR(KX) = DX_ERR DX_X_W(KX) = 1.0D0/DX_ERR DY_X(KX) = DY DY_X_ERR(KX) = DY_ERR DY_X_W(KX) = 1.0D0/DY_ERR IF ( DIST_X_ERR(KX) > 0.4 ) KX = KX - 1 END IF IF ( BAND_ARR(KP) == 'K' ) THEN KK = KK + 1 TIME_K(KK) = 1990 + (MJD - MJD_0)/365.25D0 DIST_K(KK) = DIST(KP) DIST_K_ERR(KK) = DIST_ERR(KP) DIST_K_W(KK) = 1.D0/DIST_ERR(KP) END IF IF ( J1 > 1 ) THEN IF ( BAND_ARR(KP) == 'X' .AND. BAND_ARR(KP-1) == 'S' ) THEN K_SX = K_SX + 1 TIME_SX(K_SX) = 1990 + (MJD - MJD_0)/365.25D0 DIFF_SX(K_SX) = DIST(KP) - DIST(KP-1) DIST_SX_BAND(K_SX,1) = DIST(KP) DIST_SX_BAND(K_SX,2) = DIST(KP-1) DIST_SX_BAND_ERR(K_SX,1) = DIST_ERR(KP) DIST_SX_BAND_ERR(K_SX,2) = DIST_ERR(KP-1) DIFF_SX_ERR(K_SX) = DSQRT ( DIST_ERR(KP)**2 + DIST_ERR(KP-1)**2 ) WRITE ( 6, 110 ) K_SX, 'SX', TIME_SX(K_SX), DIFF_SX(K_SX), DIFF_SX_ERR(K_SX) 110 FORMAT ( 'K= ', I2, 1X, A, ' Time: ', F7.2, ' Diff_SX: ', F5.2, ' -+ ', F5.3 ) END IF IF ( BAND_ARR(KP) == 'X' .AND. BAND_ARR(KP-1) == 'C' ) THEN K_CX = K_CX + 1 TIME_CX(K_CX) = 1990 + (MJD - MJD_0)/365.25D0 DIFF_CX(K_CX) = DIST(KP) - DIST(KP-1) DIFF_CX_ERR(K_CX) = DSQRT ( DIST_ERR(KP)**2 + DIST_ERR(KP-1)**2 ) WRITE ( 6, 110 ) K_CX, 'CX', TIME_CX(K_CX), DIFF_CX(K_CX), DIFF_CX_ERR(K_CX) END IF END IF 410 CONTINUE IUER = -1 IF ( KP > 1110 ) THEN CALL DIAGI_1E ( KP, FREQ, DIST, DIST_ERR, IUER ) END IF ! NB = LEN(BANDS) DO 420 J2=1,NB DIST_BAND(J2) = 0.0 WEI_SUM = 0.0D0 NP = 0 DO 430 J3=1,KP IF ( BAND_ARR(J3) == BANDS(J2:J2) ) THEN WEI = 1.0D0/DSQRT( DIST_ERR(J3)**2 + FLOOR_DIST_ERR**2 ) FREQ_BAND(J2) = FREQ(J3) LOG_FREQ(J2) = DLOG10 ( FREQ(J3) ) DIST_BAND(J2) = DIST_BAND(J2) + DIST(J3)*WEI WEI_SUM = WEI_SUM + WEI NP = NP + 1 END IF !! write ( 6, * ) 'j2 ', int2(j2), ' np = ', np, 'j3= ', j3, ' band_arr= ', band_arr(j3) ; call flush ( 6 ) ! %%%%%%%%%%% 430 CONTINUE DIST_BAND(J2) = DIST_BAND(J2)/WEI_SUM LOG_DIST(J2) = LOG ( DIST_BAND(J2) ) DIST_BAND_ERR(J2) = DSQRT ( 1.D0/NP/WEI_SUM**2 + FLOOR_DIST_ERR**2 ) 420 CONTINUE ! CALL REGR8 ( NB, LOG_FREQ(1), LOG_DIST(1), DR, SH, IUER ) ! WRITE ( 6, * ) 'DR= ', DR ! IV = 1 !! CALL RGRW8 ( KK, TIME_K, DIST_K, DIST_K_W, IV, MEAN_T_K, DR_K, SH_K, DR_K_ERR, SH_ERR, IUER ) CALL RGRW8 ( KX, TIME_X, DX_X, DX_X_W, IV, MEAN_X_T, DR_X_X, SH_X_X, DR_X_X_ERR, SH_X_X_ERR, IUER ) CALL RGRW8 ( KX, TIME_X, DY_X, DY_X_W, IV, MEAN_Y_T, DR_Y_X, SH_Y_X, DR_Y_X_ERR, SH_Y_X_ERR, IUER ) WRITE ( 6, 122 ) ATAN_CS ( -SH_Y_X, -SH_X_X )/DEG__TO__RAD, ( DABS(SH_X_X)*SH_Y_X_ERR + DABS(SH_Y_X)*SH_X_X_ERR )/ & & DSQRT(SH_X_X**2 + SH_Y_X**2)/DEG__TO__RAD 122 FORMAT ( 'X-band positiun angle: ', F8.1, ' -+ ', F8.1, ' deg' ) CALL RGRW8 ( KX, TIME_X, DIST_X, DIST_X_W, IV, MEAN_T_X, DR_X, SH_X, DR_X_ERR, SH_ERR, IUER ) CALL RGRW8 ( KC, TIME_C, DIST_C, DIST_C_W, IV, MEAN_T_C, DR_C, SH_C, DR_C_ERR, SH_ERR, IUER ) CALL RGRW8 ( KS, TIME_S, DIST_S, DIST_S_W, IV, MEAN_T_S, DR_S, SH_S, DR_S_ERR, SH_ERR, IUER ) WRITE ( 6, 120 ) 'K', (DIST_K(2) - DIST_K(1))/(TIME_K(2) - TIME_K(1)), (DIST_K_ERR(1) + DIST_K_ERR(2))/(TIME_K(2) - TIME_K(1)) WRITE ( 6, 120 ) 'dX', DR_X_X, DR_X_X_ERR WRITE ( 6, 120 ) 'dY', DR_Y_X, DR_Y_X_ERR WRITE ( 6, 120 ) 'Dv', DSQRT(DR_X_X**2 + DR_Y_X**2), (DABS(DR_X_X)*DR_X_X_ERR + DABS(DR_Y_X)*DR_Y_X_ERR)/ & & DSQRT(DR_X_X**2 + DR_Y_X**2) WRITE ( 6, 120 ) 'X', DR_X, DR_X_ERR WRITE ( 6, 120 ) 'C', DR_C, DR_C_ERR WRITE ( 6, 120 ) 'S', DR_S, DR_S_ERR ! DIR_ANG = ATAN_CS ( DR_Y_X, DR_X_X ) - PI__NUM DIR_ANG_ERR = (ABS(DR_X_X)*DR_Y_X_ERR + ABS(DR_Y_X)*DR_X_X_ERR)/(DR_X_X**2 + DR_Y_X**2) WRITE ( 6, 125 ) DIR_ANG/DEG__TO__RAD, DIR_ANG_ERR/DEG__TO__RAD 125 FORMAT ( 'X-band motion direction: ', F7.1, ' -+ ', F6.1, ' deg' ) ! WRITE ( 6, 130 ) 'K', 3.261*7.78*(DIST_K(2) - DIST_K(1))/(TIME_K(2) - TIME_K(1)), & & 3.261*7.78*(DIST_K_ERR(1) + DIST_K_ERR(2))/(TIME_K(2) - TIME_K(1)) WRITE ( 6, 130 ) 'X', 3.261*7.78*DR_X, 3.261*7.78*DR_X_ERR WRITE ( 6, 130 ) 'C', 3.261*7.78*DR_C, 3.261*7.78*DR_C_ERR WRITE ( 6, 130 ) 'S', 3.261*7.78*DR_S, 3.261*7.78*DR_S_ERR WRITE ( 6, 130 ) 'Dv', 3.261*7.78*DSQRT(DR_X_X**2 + DR_Y_X**2), & & 3.261*7.78*(DABS(DR_X_X)*DR_X_X_ERR + DABS(DR_Y_X)*DR_Y_X_ERR)/ & & DSQRT(DR_X_X**2 + DR_Y_X**2) ! 120 FORMAT ( A,'-band motion: ', F7.4, ' -+ ', F6.4, ' mas/yr ' ) 130 FORMAT ( A,'-band motion: ', F5.2, ' -+ ', F5.2, ' c' ) ! CALL DIAGI_1E ( KX, TIME_X, TRA_X, DIST_X_ERR, IUER ) CALL DIAGI_1E ( KX, TIME_X, DIST_X, DIST_X_ERR, IUER ) FREQ_BAND = FREQ_BAND*1.D-9 IUER = -1 CALL CORESHIFT_MOD ( NB, FREQ_BAND, DIST_BAND, DIST_BAND_ERR, & & MP, FREQ_MOD, DIST_MOD, DIST_MOD_ERR, IUER ) ! CALL DIAGI_SETDEF ( IUER, 'DIAGI_ILST', 2 ) CALL DIAGI_SETDEF ( IUER, 'DIAGI_IBST', 2 ) CALL DIAGI_SETDEF ( IUER, 'DIAGI_IPST', 5 ) CALL DIAGI_2E ( NB, FREQ_BAND, DIST_BAND, DIST_BAND_ERR, & & MP, FREQ_MOD, DIST_MOD, DIST_MOD_ERR, IUER ) CALL DIAGI_2E ( K_SX, TIME_SX, DIST_SX_BAND(1,1), DIST_SX_BAND_ERR(1,1), & & K_SX, TIME_SX, DIST_SX_BAND(1,2), DIST_SX_BAND_ERR(1,2), IUER ) DO 440 J4=1,KX TREND_X(J4) = SH_X + DR_X*(TIME_X(J4) - MEAN_T_X) TREND_X_ERR(J4) = 0.001 440 CONTINUE ! CALL DIAGI_1E ( KX, TIME_X, DIST_X, DIST_X_ERR, IUER ) CALL DIAGI_2E ( KX, TIME_X, DIST_X, DIST_X_ERR, KX, TIME_X, TREND_X, TREND_X_ERR, IUER ) ! CALL DIAGI_1E ( KC, TIME_C, DIST_C, DIST_C_ERR, IUER ) ! ! CALL DIAGI_1E ( KS, TIME_S, DIST_S, DIST_S_ERR, IUER ) ! CALL DIAGI_1E ( KC, TIME_C, DIST_C, DIST_C_ERR, IUER ) ! CALL DIAGI_1E ( NB, FREQ_BAND, DIST_BAND, DIST_BAND_ERR, IUER ) ! CALL DIAGI_1 ( NB, LOG_FREQ, LOG_DIST, IUER ) ! CALL DIAGI_1E ( K_SX, TIME_SX, DIFF_SX, DIFF_SX_ERR, IUER ) ! CALL DIAGI_1E ( K_CX, TIME_CX, DIFF_CX, DIFF_CX_ERR, IUER ) ! END PROGRAM CORE_SHIFT !#! ! ! ------------------------------------------------------------------------ ! SUBROUTINE CORESHIFT_MOD ( NB, FREQ, DIST_OBS, DIST_ERR, & & MP, FREQ_MOD, DIST_MOD, DIST_MOD_ERR, IUER ) IMPLICIT NONE INTEGER*4 NB, MP, IUER REAL*8 FREQ(NB), DIST_OBS(NB), DIST_ERR(NB), & & FREQ_MOD(MP), DIST_MOD(MP), DIST_MOD_ERR(MP) REAL*8 FREQ_PLT INTEGER*4 MT PARAMETER ( MT = 128 ) REAL*8 EXP_MIN, EXP_MAX, EXP_VAL, EXP_BEST, & & SHF_MIN, SHF_MAX, SHF_VAL, SHF_BEST, SH, DR, RES, WW, & & WRMS_SQ, WRMS_SQ_MIN REAL*8 DIST_OBS_TRY(NB), EXP_TRY(MT), SHF_TRY(MT), EXP_AVR, EXP_SQR, & & EXP_RMS REAL*8, EXTERNAL :: RGAUSS INTEGER*4 J1, J2, J3, J4, J5, J6, IS, NF, NEX, NSH ! FREQ_PLT = 50.0 EXP_MIN = -1.5 EXP_MAX = -0.5 SHF_MIN = -6.0 SHF_MAX = -2.0 EXP_AVR = 0.0D0 EXP_SQR = 0.0D0 NEX = 512 NSH = 512 NF = 1 IS = 8105156 DO 410 J1=1,MT WW = 0.0D0 WRMS_SQ_MIN = 1.0D8 DO 420 J2=NF,NB WW = WW + 1.0D0/DIST_ERR(J2) DIST_OBS_TRY(J2) = DIST_OBS(J2) + RGAUSS ( IS, DIST_ERR(J2) ) IF ( J1 == MT ) THEN DIST_OBS_TRY(J2) = DIST_OBS(J2) END IF 420 CONTINUE DO 430 J3=1,NEX EXP_VAL = EXP_MIN + (J3-1)*(EXP_MAX - EXP_MIN)/(NEX-1) DO 440 J4=1,NSH SHF_VAL = SHF_MIN + (J4-1)*(SHF_MAX - SHF_MIN)/(NSH-1) RES = 0.0D0 DO 450 J5=NF,NB DIST_MOD(J5) = DIST_OBS_TRY(NB) + SHF_VAL*FREQ(J5)**EXP_VAL RES = RES + (DIST_OBS_TRY(J5) - DIST_MOD(J5))**2/DIST_ERR(J5) 450 CONTINUE WRMS_SQ = RES/WW IF ( WRMS_SQ < WRMS_SQ_MIN ) THEN WRMS_SQ_MIN = WRMS_SQ SHF_BEST = SHF_VAL EXP_BEST = EXP_VAL END IF 440 CONTINUE 430 CONTINUE EXP_TRY(J1) = EXP_BEST SHF_TRY(J1) = SHF_BEST EXP_AVR = EXP_AVR + EXP_TRY(J1) EXP_SQR = EXP_SQR + EXP_TRY(J1)**2 410 CONTINUE EXP_AVR = EXP_AVR/MT EXP_RMS = DSQRT ( EXP_SQR/MT - EXP_AVR**2 )/DSQRT(2.0D0) ! WRITE ( 6, 110 ) EXP_AVR, EXP_RMS, WRMS_SQ_MIN 110 FORMAT ( 'EXP: ', F5.2, ' -+ ', F5.2, ' wrms= ', f10.3 ) ! DO 460 J6=1,MP FREQ_MOD(J6) = FREQ(1) + (J6-1)*(FREQ_PLT - FREQ(1))/(MP-1) DIST_MOD(J6) = DIST_OBS(NB) + SHF_BEST*FREQ_MOD(J6)**EXP_BEST DIST_MOD_ERR(J6) = WRMS_SQ_MIN 460 CONTINUE ! CALL ERR_LOG ( 0, IUER ) RETURN END SUBROUTINE CORESHIFT_MOD !#!