* Last processed by NICE on 12-Jun-2000 15:54:00 * Customized for : IEEE, LINUX, UNIX, MOTIF, F77 PROGRAM SHIFT_AND_ADD * ========================= * * This is the main program of the gildas task MAKPLAN. * It acceeds to SIC variables defined by NIC , when reading a file, and * to 2 shared memory area which contain the subscan data and the record data. * Then it calls the subroutine MAKPLAN which grids data. * INCLUDE 'inc:memory.inc' INCLUDE 'inc:format.inc' INCLUDE 'inc:xpar.inc' INCLUDE 'const.inc' INCLUDE 'parameter.inc' INCLUDE 'makplan.inc' INCLUDE 'radec.inc' * * real*8 * REAL*8 W_MIN,W_MAX REAL*8 ALIM,ALOW REAL*8 GONS(100,4), BOUND(5) REAL*8 DRA,DDEC,RA_FROM,DEC_FROM,RA_DAY,DEC_DAY,SITE_LATITUDE REAL*8 DAZ,DEL,DLST C Shared memory parameters INTEGER*4 POINTER INTEGER*4 IADDRRECORD ! address of shared memory area INTEGER*4 IP_AZOFF,IP_ELOFF,IP_RECEIVER_FLAG INTEGER*4 IPSUBSCAN,IP_EPOCH,IPSUBREC INTEGER*4 IPRECORD,IP_LST,IP_EL,IP_SIGNAL,IP_SCAN_COORD INTEGER*4 IP,IPMAP,IP_PARAM INTEGER*4 IPMAPD_AZ,IPMAPW_AZ,IPMAPD_EL,IPMAPW_EL INTEGER*4 IP_NEWSIGNAL, IP_NEWXCOORD INTEGER*4 IP_LST_MAP,IP_LST_ROW INTEGER*4 IP_EL_MAP,IP_EL_ROW INTEGER*4 IP_WEIGHT_MAP,IP_WEIGHT_ROW INTEGER*4 IP_OFFX_MAP,IP_OFFY_MAP,IP_ANGLE_MAP INTEGER*4 IADDR_MAP, IADDR INTEGER*4 IP_SUPPORT, IP_SUPPORT_CHAN INTEGER SIC_GETVM INTEGER DATA_SHMID ! shared memory identifier INTEGER MAX_SEG , NSUP PARAMETER (MAX_SEG=100) * * INTEGER IER ! error status INTEGER NMAP REAL FREQ LOGICAL RECEIVER_FLAG INTEGER I,J,IJ,II,ICHAN,IMAP,IFIRST,ILAST,ISUBSCAN INTEGER NP_SUBSCAN REAL VAL,WEI REAL MAXMIN(8) REAL UNBAL,WOBBLER_THROW LOGICAL ERROR REAL P_ANGLE, CONV_ANGLE REAL OFFX,OFFY LOGICAL SUPPORT_OPT INTEGER REF_CHAN INTEGER N_RECEIVERS,NPARAM_RECEIVER,NPARAM_SUBSCAN INTEGER FIELD_ROTATION * INTEGER SCAN_TYPE,IFUN INTEGER BLC(4), TRC(4) * REAL OLAM,OBET COMMON /SHIFT_AND_ADD_COMMON/ $DRA,DDEC,RA_FROM,DEC_FROM,RA_DAY,DEC_DAY,SITE_LATITUDE,DLST, $OLAM,OBET SAVE /SHIFT_AND_ADD_COMMON/ * DATA BLC /0,0,0,0/, TRC /0,0,0,0/ * * * Read gildas parameters * CALL GILDAS_OPEN CALL GILDAS_INTE('DATA_SHMID$',DATA_SHMID,1) * CALL GILDAS_INTE('NSUBSCAN$',NSUBSCAN,1) CALL GILDAS_INTE('NRECORD$',NRECORD,1) CALL GILDAS_INTE('SCAN_TYPE$',SCAN_TYPE,1) * IF (DATA_SHMID.EQ.0 .OR. NSUBSCAN.EQ.0 .OR. NRECORD.EQ.0) THEN WRITE(6,'(/a,/)') $ 'No valid data in memory. Cannot run SAA' CALL GILDAS_CLOSE STOP ENDIF * IF (SCAN_TYPE.NE.2) THEN WRITE(6,'(/a,/)') $ 'Observing mode is not "ON_THE_FLY". Cannot run SAA.' CALL GILDAS_CLOSE STOP ENDIF * CALL GILDAS_INTE('SCAN_NUMBER$',SCAN_NUMBER,1) CALL GILDAS_INTE('NCHAN$',NCHANNEL,1) CALL GILDAS_LOGI('JANSKY_FLAG$',JANSKY_FLAG,1) CALL GILDAS_REAL('FREQUENCY$',FREQ,1) CALL GILDAS_CHAR('SNAM$',SNAM) CALL GILDAS_REAL('SLAM$',SLAM,1) CALL GILDAS_REAL('SBET$',SBET,1) CALL GILDAS_DBLE('EPOCH$',EPOCHFR,1) CALL GILDAS_DBLE('SITE_LATITUDE$',SITE_LATITUDE,1) CALL GILDAS_REAL('OAZM$',OAZM,1) CALL GILDAS_REAL('OELV$',OELV,1) CALL GILDAS_REAL('XINC$',XINC,1) CALL GILDAS_REAL('YINC$',YINC,1) CALL GILDAS_REAL('SINT$',SINT,1) IF (SINT.EQ.0) THEN WRITE(6,'(/a,/)') $ 'Observing mode is not ON_THE_FLY. Cannot run SAA.' CALL GILDAS_CLOSE STOP ENDIF CALL GILDAS_REAL('ZIGZAG$',ZIGZAG,1) CALL GILDAS_LOGI('FLIP$',FLIP,1) CALL GILDAS_LOGI('SMOOTH_EL$',SMOOTH_EL,1) CALL GILDAS_INTE('FUNCTION$',IFUN,1) CALL GILDAS_REAL('BEAMSEP$',WOBBLER_THROW,1) CALL GILDAS_REAL('UNBAL$',UNBAL,1) CALL GILDAS_REAL('PA$',P_ANGLE,1) CALL GILDAS_REAL('OFFX$',OFFX,1) CALL GILDAS_REAL('OFFY$',OFFY,1) CALL GILDAS_LOGI('SAA_MASK$',SUPPORT_OPT,1) CALL GILDAS_CHAR('TELESCOPE$',TELESCOPE) CALL GILDAS_INTE('REF_CHAN$',REF_CHAN,1) CALL GILDAS_INTE('FIELD_ROTATION$',FIELD_ROTATION,1) * CALL GILDAS_CLOSE * IF (WOBBLER_THROW.LE.1) THEN WRITE (6,'(/a,/)') 'No valid BEAMSEP. Cannot run SAA.' STOP ENDIF * * Conversion FREQUENCY = FREQ*1.D9 ! frequency in Herz * * * * Link to data shared memory C+WIN32 c CALL OPEN_SHARED_MEMORY(DATA_SHMID) C-WIN32 CALL ATTACH_MEMORY(DATA_SHMID,IADDRRECORD) IP = POINTER(IADDRRECORD,MEMORY) CALL R4TOR4(MEMORY(IP),N_RECEIVERS,1) CALL R4TOR4(MEMORY(IP+1),NPARAM_RECEIVER,1) CALL R4TOR4(MEMORY(IP+2),NPARAM_SUBSCAN,1) IP_PARAM = IP+3 IP_AZOFF = IP_PARAM IP_ELOFF = IP_PARAM + NCHANNEL IP_RECEIVER_FLAG = IP_PARAM + 2*NCHANNEL IPSUBSCAN = IP_PARAM + NPARAM_RECEIVER*NCHANNEL IP_EPOCH = IPSUBSCAN + NSUBSCAN IPSUBREC = IPSUBSCAN + 3*NSUBSCAN IPRECORD = IPSUBSCAN + 4*NSUBSCAN IP_LST = IPRECORD IP_EL = IPRECORD + 2*NRECORD IP_SIGNAL = IPRECORD + 3*NRECORD IP_SCAN_COORD = IP_SIGNAL + N_RECEIVERS*NRECORD IP_SUPPORT = IP_SCAN_COORD + 2*NRECORD + $4*NCHANNEL*MAX_SEG+1 ! ip_support points on ekh/saa support * * Zigzag correction IF (ZIGZAG.NE.0) CALL ZIGZAG_CORRECTION( $MEMORY(IP_SCAN_COORD),MEMORY(IPSUBREC)) * * Allocate memory for regridded data in azimuth and in elevation (data and * weight) : mapd_a(nlin,ncol,nchannel),mapw_a(),mapw_e() * First find dimension of map (according to xcoord, ycoord and map_cell) CALL MAXMIN_COORD ( $MEMORY(IP_SCAN_COORD),NRECORD,XMIN,XMAX,MEMORY(IP_SIGNAL)) IF (XMIN-OAZM*3600.LT.SINT) XMIN = OAZM*3600. NCOL = (XMAX-XMIN)/XINC +1 * NROW = NSUBSCAN CALL MAXMIN_COORD( $MEMORY(IP_SCAN_COORD+NRECORD),NRECORD,YMIN,YMAX, $MEMORY(IP_SIGNAL)) IF (YINC.NE.SINT) THEN NROW = (YMAX-YMIN)/YINC +1 ENDIF NMAP = NCOL*NROW*NCHANNEL * IER = SIC_GETVM(2*NMAP+5*NROW*NCOL+4*NCOL,IADDR_MAP) IF (IER.NE.1) THEN WRITE(6,*)' *** Unable to allocate memory (1)***' STOP ENDIF IPMAP = POINTER(IADDR_MAP,MEMORY) IPMAPD_AZ = IPMAP ! length = NCOL IPMAPW_AZ = IPMAPD_AZ + NCOL ! length = NCOL IPMAPD_EL = IPMAPW_AZ + NCOL ! length = NCOL*NROW*NCHANNEL IPMAPW_EL = IPMAPD_EL + NMAP ! length = NCOL*NROW*NCHANNEL IP_LST_MAP = IPMAPW_EL + NMAP ! length = NCOL*NROW IP_LST_ROW = IP_LST_MAP + NCOL*NROW ! length = NCOL IP_EL_MAP = IP_LST_ROW + NCOL ! length = NCOL*NROW IP_EL_ROW = IP_EL_MAP + NCOL*NROW ! length = NCOL IP_OFFX_MAP = IP_EL_ROW + NCOL ! length = NCOL*NROW IP_OFFY_MAP = IP_OFFX_MAP + NCOL*NROW ! length = NCOL*NROW IP_ANGLE_MAP = IP_OFFY_MAP + NCOL*NROW ! length = NCOL*NROW IP_WEIGHT_MAP= IPMAPW_EL IP_WEIGHT_ROW= IPMAPW_AZ * * * -- Initialize convolving function * CALL FUNGEN(FUNC,IFUN) * * -- Test if support exists CALL R4TOR4(MEMORY(IP_SUPPORT),NSUP,1) IF (SUPPORT_OPT) THEN IF (NSUP.EQ.0) THEN WRITE(6,'(/,a)') $ '*** No support available ***' SUPPORT_OPT = .FALSE. ENDIF ENDIF * * -- Initialize LST and EL map * CALL LST_EL_MAP(MEMORY(IP_LST),MEMORY(IP_LST_MAP), $MEMORY(IP_LST_ROW),MEMORY(IP_EL),MEMORY(IP_EL_MAP), $MEMORY(IP_EL_ROW),MEMORY(IP_WEIGHT_MAP),MEMORY(IP_WEIGHT_ROW), $MEMORY(IP_SCAN_COORD),MEMORY(IP_SCAN_COORD+NRECORD), $MEMORY(IPSUBREC),IFUN) * * Compute RA_DAY and DEC_DAY CALL R4TOR8(MEMORY(IP_EPOCH),EPOCHTO,1) CALL R4TOR8(MEMORY(IP_LST+NRECORD/2),DLST,1) RA_FROM = SLAM + OLAM DEC_FROM = SBET + OBET DRA = 0.D0 DDEC = 0.D0 CALL SETUP CALL PREC(3,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA_FROM,DEC_FROM) * SITE_LATITUDE = SITE_LATITUDE /DEG_TO_RAD * SINDEC = DSIN(DEC_DAY*DEG_TO_RAD) ! init for radec COSDEC = DCOS(DEC_DAY*DEG_TO_RAD) ! init for radec SINOL = DSIN(SITE_LATITUDE*DEG_TO_RAD) ! init for radec COSOL = DCOS(SITE_LATITUDE*DEG_TO_RAD) ! init for radec * * * -- Initialize OFFX and OFFY map * DRA = OFFX DDEC = OFFY DO J=1,NROW DO I=1,NCOL IJ = I + (J-1)*NCOL - 1 CALL R4TOR8(MEMORY(IP_LST_MAP+IJ),DLST,1) CALL RADEC_TO_AZEL(DRA,DDEC,DLST,SITE_LATITUDE, $ RA_DAY, DEC_DAY,DAZ,DEL) CALL R8TOR4(DAZ,MEMORY(IP_OFFX_MAP+IJ),1) CALL R8TOR4(DEL,MEMORY(IP_OFFY_MAP+IJ),1) ENDDO ENDDO * PA * -- Initialize angle map * DRA = SIN(P_ANGLE*DEG_TO_RAD)*WOBBLER_THROW DDEC = COS(P_ANGLE*DEG_TO_RAD)*WOBBLER_THROW DO J=1,NROW DO I=1,NCOL IJ = I + (J-1)*NCOL - 1 CALL R4TOR8(MEMORY(IP_LST_MAP+IJ),DLST,1) CALL RADEC_TO_AZEL(DRA,DDEC,DLST,SITE_LATITUDE, $ RA_DAY, DEC_DAY,DAZ,DEL) CONV_ANGLE = ATAN2(DEL,DAZ)/DEG_TO_RAD CALL R4TOR4(CONV_ANGLE,MEMORY(IP_ANGLE_MAP+IJ),1) ENDDO ENDDO C C C -- Raz array DO I = 0,NCOL*NROW*NCHANNEL-1 CALL R4TOR4(0,MEMORY(IPMAPD_EL + I),1) CALL R4TOR4(0,MEMORY(IPMAPW_EL + I),1) ENDDO * * * C CALCULATE NCHANNEL MAPS C ------------------------ * DO ICHAN=1,NCHANNEL IMAP = (ICHAN-1)*NCOL*NROW ILAST = 0 CALL R4TOR4(MEMORY(IP_RECEIVER_FLAG+ICHAN-1), $ RECEIVER_FLAG,1) IF (.NOT.RECEIVER_FLAG) THEN DO ISUBSCAN=1,NSUBSCAN CALL R4TOR4(MEMORY(IPSUBREC+ISUBSCAN-1),NP_SUBSCAN,1) IFIRST = ILAST +1 ILAST = IFIRST + NP_SUBSCAN - 1 IF (NP_SUBSCAN.NE.0) THEN IER = SIC_GETVM(4*NP_SUBSCAN,IADDR) IF (IER.NE.1) THEN WRITE(6,*) $ ' *** Unable to allocate memory (2)***' STOP ENDIF IP_NEWSIGNAL = POINTER(IADDR,MEMORY) IP_NEWXCOORD = IP_NEWSIGNAL + 2*NP_SUBSCAN CALL SHIFT_ADD ( $ MEMORY(IP_SIGNAL+IFIRST-1+(ICHAN-1)*NRECORD), $ MEMORY(IP_SCAN_COORD+IFIRST-1), $ MEMORY(IP_NEWSIGNAL),MEMORY(IP_NEWXCOORD), $ NP_SUBSCAN,WOBBLER_THROW,UNBAL) * c -- c Nasmith ---> Az/El rotation c Output XGRID,YGRID (rotated XGRID,YGRID) and ZGRID (ori. XGRID) c -- * * ! The subroutine azcon regrids the data in azimuth, and * ! places the resultant scan into the array row_d. * * CALL AZCON (MEMORY(IPMAPD_AZ), $ MEMORY(IPMAPW_AZ), $ MEMORY(IP_NEWSIGNAL),MEMORY(IP_NEWXCOORD), $ 2*NP_SUBSCAN,IFUN) * * * * ! The subroutine ELCON makes the az-el map from the azimuth * ! scans by convolution of the data onto the number grid of * ! the map array 'mapd_e' * ! In the process it noise-filters in * ! the elevation direction and regrids the data to allow * ! for an elevation offset. * ! Input : * ! - mapd_a, mapw_a: data and weight after gridding in azimuth * ! Output : * ! - mapd_e, mapw_e: data and weight after gridding in elevation * IF (SMOOTH_EL.OR.(NROW.NE.NSUBSCAN)) THEN CALL ELCON (MEMORY(IPMAPD_AZ), $ MEMORY(IPMAPW_AZ), $ MEMORY(IPMAPD_EL+IMAP), $ MEMORY(IPMAPW_EL+IMAP), $ MEMORY(IP_SCAN_COORD+NRECORD+IFIRST-1),IFUN) ELSE * !copy row_d and row_w in maps_el and mapw (nrow=nsubscan) II = (ISUBSCAN-1)*NCOL + IMAP DO I=0,NCOL-1 * MAP_D(II+I) = ROW_D(I) * ROW_W(I) * MAP_W(II+I) = ROW_W(I) CALL R4TOR4(MEMORY(IPMAPD_AZ+I),VAL,1) CALL R4TOR4(MEMORY(IPMAPW_AZ+I),WEI,1) VAL = VAL*WEI CALL R4TOR4(VAL,MEMORY(IPMAPD_EL+II+I),1) CALL R4TOR4(WEI,MEMORY(IPMAPW_EL+II+I),1) ENDDO ENDIF * * CALL FREE_VM(4*NP_SUBSCAN,IADDR) ENDIF ENDDO * * --- Mask data * IF (SUPPORT_OPT) THEN IP_SUPPORT_CHAN = IP_SUPPORT+1 + 2*NSUP*(ICHAN-1) CALL INIT_SUPPORT_MASK (GONS,BOUND,MAX_SEG,NSUP, $ MEMORY(IP_SUPPORT_CHAN),MEMORY(IP_SUPPORT_CHAN+NSUP)) CALL SUPPORT_MASK(MEMORY(IPMAPD_EL+IMAP), $ XMIN,XINC,YMIN,YINC, $ NCOL,NROW,GONS,NSUP,BOUND, $ WOBBLER_THROW) ELSE CALL ANGLE_MASK(MEMORY(IPMAPD_EL+IMAP), $ MEMORY(IP_AZOFF),MEMORY(IP_ELOFF), $ XMIN,XINC,YMIN,YINC, $ MEMORY(IP_EL_MAP), $ MEMORY(IP_OFFX_MAP),MEMORY(IP_OFFY_MAP), $ MEMORY(IP_ANGLE_MAP), $ NCOL,NROW,ICHAN,REF_CHAN,WOBBLER_THROW,TELESCOPE, $ FIELD_ROTATION) ENDIF * * C --- Compute maximum and minimum weight W_MIN = 0 W_MAX = 0 CALL WMAX(MEMORY(IPMAPD_EL+IMAP), $ MEMORY(IPMAPW_EL+IMAP), $ W_MAX,W_MIN,BLANKING,NCOL*NROW) * C --- Normalize : C divide by weight C flagge data if weight < 0.1*w_max (original value: 0.2*w_max) C or weight < 0.5*w_min (weight negative) C ALIM = 0.1*W_MAX ALOW = 0.5*W_MIN CALL MNORM(MEMORY(IPMAPD_EL+IMAP), $ MEMORY(IPMAPW_EL+IMAP), $ ALIM,ALOW,BLANKING,NCOL*NROW,1.) * ELSE ! flagged channel DO I = 0,NCOL*NROW-1 CALL R4TOR4(BLANKING,MEMORY(IPMAPD_EL+IMAP+I),1) ENDDO ENDIF ENDDO * * restore xcoord (remove zigzag correction) IF (ZIGZAG.NE.0) CALL REMOVE_ZIGZAG_CORRECTION( $MEMORY(IP_SCAN_COORD),MEMORY(IPSUBREC)) * * C --- make GILDAS image file * *RN-17SEP96 CALL MAXMIN_IMAGE (MEMORY(IPMAPD_EL+IMAP),NCOL,NROW, CALL MAXMIN_IMAGE (MEMORY(IPMAPD_EL),NCOL,NROW, $NCHANNEL,MAXMIN) X_NDIM = 3 IF (.NOT.JANSKY_FLAG) THEN X_UNIT = ' Counts' ELSE IF (MAXMIN(2).LT.2) THEN X_UNIT = ' mJy' DO ICHAN = 1,NCHANNEL IMAP = (ICHAN-1)*NCOL*NROW CALL NORM3 (MEMORY(IPMAPD_EL+IMAP), $ MEMORY(IPMAPW_EL+IMAP), $ ALIM,ALOW,BLANKING,NCOL*NROW) ENDDO *RN-17SEP96 CALL MAXMIN_IMAGE (MEMORY(IPMAPD_EL+IMAP), CALL MAXMIN_IMAGE (MEMORY(IPMAPD_EL), $ NCOL,NROW,NCHANNEL,MAXMIN) ELSE X_UNIT = ' Jy' ENDIF ENDIF X_DIM(1) = NCOL X_DIM(2) = NROW X_DIM(3) = NCHANNEL X_DIM(4) = 1 X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4) * X_REF1 = 1 X_VAL1 = -ABS(XMIN)*SEC_TO_RAD X_INC1 = XINC*SEC_TO_RAD * X_REF2 = 1 IF (SMOOTH_EL.OR.(NROW.NE.NSUBSCAN)) THEN X_VAL2 = -ABS(OELV)*DEG_TO_RAD X_INC2 = YINC*SEC_TO_RAD ELSE X_VAL2 = OELV*DEG_TO_RAD X_INC2 = -SIGN(1.,OELV)*YINC*SEC_TO_RAD ENDIF * X_REF3 = 0 X_VAL3 = 0 X_INC3 = 1 * X_REF4 = 0 X_VAL4 = 0 X_INC4 = 1 * X_TYPE = 'GILDAS_IMAGE' X_NAME = SNAM X_FORM = FMT_R4 X_GENE = 29 X_SYST = 'AZIMUTHAL' X_BLAN = 2 X_BVAL = BLANKING X_EVAL = 0. X_EXTR = 10 X_RMIN = MAXMIN(1) X_RMAX = MAXMIN(2) X_MIN1 = MAXMIN(3) X_MIN2 = MAXMIN(4) X_MIN3 = MAXMIN(5) X_MAX1 = MAXMIN(6) X_MAX2 = MAXMIN(7) X_MAX3 = MAXMIN(8) X_MIN4 = 1 X_MAX4 = 1 * X_DESC = 18 X_POSI = 12 X_RA = SLAM*DEG_TO_RAD X_DEC = SBET*DEG_TO_RAD X_PROJ = 9 X_SPEC = 12 X_RESO = 3 X_FREQ = FREQUENCY/1.D6 ! in Mhz * SNAM = ' ' WRITE (SNAM(1:4),'(i4.4)') SCAN_NUMBER CALL SIC_LOWER (SNAM) CALL SIC_PARSEF(SNAM,X_FILE,' ','.saa') CALL GDF_GEIS (X_ISLO,ERROR) ! reserve slot CALL GDF_WRITX (X_ISLO,ERROR) CALL GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR) CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR) IP = POINTER (X_ADDR,MEMORY) CALL R4TOR4(MEMORY(IPMAPD_EL),MEMORY(IP), $NCOL*NROW*NCHANNEL) CALL GDF_FRIS (X_ISLO,ERROR) CALL FREE_VM(2*NMAP+2*NCOL,IADDR_MAP) * 99 CONTINUE STOP END * *---------------------------------------------------------------- * * SHIFT_ADD() * *---------------------------------------------------------------- * C+HPUX c$OPTIMIZE LEVEL1 C-HPUX SUBROUTINE SHIFT_ADD (SIGNAL,XCOORD,NEW_SIGNAL,NEW_XCOORD, $NP_SUBSCAN,WOBBLER_THROW,UNBAL) * INCLUDE 'parameter.inc' * REAL SIGNAL(1),XCOORD(1),NEW_SIGNAL(1),NEW_XCOORD(1) REAL WOBBLER_THROW INTEGER NP_SUBSCAN REAL UNBAL * INTEGER I * DO I=1,NP_SUBSCAN IF (ABS(SIGNAL(I)-BLANKING).GT.0.5) THEN NEW_SIGNAL(I) = SIGNAL(I)*UNBAL NEW_SIGNAL(I+NP_SUBSCAN) = -SIGNAL(I)*UNBAL ELSE NEW_SIGNAL(I) = BLANKING NEW_SIGNAL(I+NP_SUBSCAN) = BLANKING ENDIF NEW_XCOORD(I) = XCOORD(I) + WOBBLER_THROW*10/2. NEW_XCOORD(I+NP_SUBSCAN) = XCOORD(I)-WOBBLER_THROW*10/2. ENDDO RETURN END *---------------------------------------------------------------- * * ANGLE_MASK() * *---------------------------------------------------------------- SUBROUTINE ANGLE_MASK (MAPD,OA,OE,XMIN,XINC,YMIN,YINC,EL, $OFFX,OFFY,ANGLE,NCOL,NROW,ICHAN,REF_CHAN, $WOBBLER_THROW,TELESCOPE,FIELD_ROTATION) * INCLUDE 'const.inc' INCLUDE 'parameter.inc' * INTEGER NCOL,NROW,ICHAN,REF_CHAN REAL MAPD(1),OA(1),OE(1) REAL EL(NCOL,NROW) REAL OFFX(NCOL,NROW),OFFY(NCOL,NROW),ANGLE(NCOL,NROW) REAL XMIN,XINC,YMIN,YINC REAL WOBBLER_THROW CHARACTER*40 TELESCOPE INTEGER FIELD_ROTATION * REAL*8 XCOO,YCOO REAL COS_EL,SIN_EL,COTG_ANGLE,DIST INTEGER I,J LOGICAL NULL_ANGLE * DO J=1,NROW DO I=1,NCOL IF (ABS(EL(I,J)-BLANKING).GT.1) THEN COS_EL = COS(FIELD_ROTATION*EL(I,J)*DEG_TO_RAD) SIN_EL = SIN(FIELD_ROTATION*EL(I,J)*DEG_TO_RAD) IF (TELESCOPE(1:3).EQ.'CSO') THEN COS_EL = 1. SIN_EL = 0 ENDIF NULL_ANGLE = .TRUE. IF (ABS(ANGLE(I,J)).GT.1-06) THEN NULL_ANGLE = .FALSE. COTG_ANGLE = COS(ANGLE(I,J)*DEG_TO_RAD) / $ SIN(ANGLE(I,J)*DEG_TO_RAD) ENDIF * YCOO = YMIN + (J-1)*YINC $ +(OA(ICHAN)-OA(REF_CHAN))*SIN_EL $ -(OE(ICHAN)-OE(REF_CHAN))*COS_EL XCOO = XMIN + (I-1)*XINC $ -(OA(ICHAN)-OA(REF_CHAN))*COS_EL $ -(OE(ICHAN)-OE(REF_CHAN))*SIN_EL IF (.NOT.NULL_ANGLE) THEN DIST = OFFX(I,J) +(YCOO-OFFY(I,J))*COTG_ANGLE -XCOO IF (ABS(DIST).GT.WOBBLER_THROW/2) $ MAPD(I+(J-1)*NCOL) = BLANKING ELSE MAPD(I+(J-1)*NCOL) = BLANKING ENDIF ELSE MAPD(I+(J-1)*NCOL) = BLANKING ENDIF ENDDO ENDDO RETURN END *---------------------------------------------------------------- * * SUPPORT_MASK() * *---------------------------------------------------------------- SUBROUTINE SUPPORT_MASK (MAPD,XMIN,XINC,YMIN,YINC, $NCOL,NROW,GONS,PNT,BOUND,WOBBLER_THROW) * INCLUDE 'const.inc' INCLUDE 'parameter.inc' * INTEGER MAX_SEG PARAMETER (MAX_SEG=100) REAL*8 GONS(MAX_SEG,4),BOUND(5) REAL MAPD(1) REAL XMIN,XINC,YMIN,YINC INTEGER NCOL,NROW,PNT REAL WOBBLER_THROW * REAL*8 XCOO,YCOO INTEGER I,J LOGICAL IN,GR8_IN INTEGER FIRST_IN,LAST_IN,N_IN * DO J=1,NROW FIRST_IN = NCOL+1 LAST_IN = 0 DO I=1,NCOL YCOO = YMIN + (J-1)*YINC XCOO = XMIN + (I-1)*XINC * Test if inside support IN = GR8_IN(XCOO,YCOO,MAX_SEG,PNT-1,GONS,BOUND) IF (IN .AND. I.LT.FIRST_IN) FIRST_IN = I IF (IN .AND. I.GT.LAST_IN) LAST_IN = I ENDDO N_IN = WOBBLER_THROW/XINC DO WHILE (LAST_IN-FIRST_IN+1 .GT. N_IN) FIRST_IN = FIRST_IN+1 IF (LAST_IN-FIRST_IN+1 .GT. N_IN) $ LAST_IN = LAST_IN-1 ENDDO DO I=1,FIRST_IN-1 MAPD(I+(J-1)*NCOL) = BLANKING ENDDO DO I=LAST_IN+1,NCOL MAPD(I+(J-1)*NCOL) = BLANKING ENDDO ENDDO END *---------------------------------------------------------------- * * INIT_SUPPORT_MASK() * *---------------------------------------------------------------- * SUBROUTINE INIT_SUPPORT_MASK (GONS,BOUND,MAX_SEG,PNT, $XU,YU) * * INTEGER PNT INTEGER MAX_SEG REAL*8 GONS(MAX_SEG,4),BOUND(5) * REAL*4 XU(1),YU(1) INTEGER I * * Initialize GONS and BOUND arrays (used by routine gr8_in) *---------------------------------------------------------- DO I=1,PNT GONS(I,1) = XU(I) GONS(I,2) = YU(I) IF (I.LT.PNT) THEN GONS(I,3) = XU(I+1)-XU(I) GONS(I,4) = YU(I+1)-YU(I) ENDIF ENDDO BOUND(2)= GONS(1,1) BOUND(3)= GONS(1,1) BOUND(4)= GONS(1,2) BOUND(5)= GONS(1,2) DO I=2,PNT BOUND(2) = MIN(BOUND(2),GONS(I,1)) ! xmin BOUND(3) = MAX(BOUND(3),GONS(I,1)) ! xmax BOUND(4) = MIN(BOUND(4),GONS(I,2)) ! ymin BOUND(5) = MAX(BOUND(5),GONS(I,2)) ! ymax ENDDO RETURN END *---------------------------------------------------------------- * * LST_EL_MAP() * *---------------------------------------------------------------- SUBROUTINE LST_EL_MAP(LST_RECORD,LST_MAP,LST_ROW, $EL_RECORD,EL_MAP,EL_ROW,WEIGHT_MAP,WEIGHT_ROW, $XSCAN_COORD,YSCAN_COORD,NSUBREC,IFUN) C C************************************************************** C This subroutine computes LST and EL corresponding to the (Az,EL) map C C C Input: C lst[NRECORD] C el[NRECORD] C xscan_coord[NRECORD] : x coordinates C yscan_coord[NRECORD] : y coordinates C np_subscan : number of records in subscan C func : convolving function ( ~ sin(x)/x ) C C Output: C lst_map : lst map C el_map : lst map C weight_map : sum of weight C C************************************************************** * INCLUDE 'const.inc' INCLUDE 'makplan.inc' INCLUDE 'parameter.inc' * * REAL*4 LST_RECORD(1), MAGI, DIAMETER REAL*4 LST_MAP(1) REAL*4 LST_ROW(1) REAL*4 EL_RECORD(1) REAL*4 EL_MAP(1) REAL*4 EL_ROW(1) REAL*4 XSCAN_COORD(1) REAL*4 YSCAN_COORD(1) REAL*4 WEIGHT_MAP(1) REAL*4 WEIGHT_ROW(1) CC REAL*4 DATAR(1) * INTEGER IFUN INTEGER NSUBREC(1) * REAL*4 FUNWIDX,TOTWIDX,FUNWIDY,TOTWIDY INTEGER IFIRST,ILAST,I,J,IFUNC,ISUBSCAN,IJ REAL*4 POSX,POSY REAL*4 X,Y,Z INTEGER NP_SUBSCAN,IFIRST_REC * DO I=1,NROW DO J=1,NCOL IJ = J+(I-1)*NCOL LST_MAP(IJ) = 0 EL_MAP(IJ) = 0 WEIGHT_MAP(IJ) = 0 ENDDO ENDDO C C C The unit length of the grid is XINC C Compute function window width in the new system unit : C funwid = lambda/(1.67*diam_telescope)/xinc C funwid = 1.34 if freq=230ghz, diam_telescope=30m and xinc=4arcsec C C * MAGI = 0 IF (IFUN.EQ.1) MAGI = 1.6806492007530D0 ! Parabolic IF (IFUN.EQ.2) MAGI = 1.6982249549743D0 ! Gaussian IF (IFUN.EQ.3) MAGI = 1.8363197305085D0 ! Linear IF (MAGI.EQ.0) MAGI = 1.6574002064571D0 ! No taper C DIAMETER = 30 IF (TELESCOPE(1:3).EQ.'CSO') DIAMETER = 10 IF (TELESCOPE(1:3).EQ.'HHT') DIAMETER = 10 FUNWIDX = VLIGHT/FREQUENCY /(MAGI*DIAMETER) $/(XINC*SEC_TO_RAD) IF ( FUNWIDX.LT.1.D0) FUNWIDX = 1.D0 TOTWIDX = 4.D0 * FUNWIDX ! total window width FUNWIDY = VLIGHT/FREQUENCY /(MAGI*DIAMETER) $/(YINC*SEC_TO_RAD) IF ( FUNWIDY.LT.1.D0) FUNWIDY = 1.D0 TOTWIDY = INT(0.9999*SINT/YINC) * * -- Loop on subscan IFIRST_REC = 1 DO ISUBSCAN=1,NSUBSCAN NP_SUBSCAN = NSUBREC(ISUBSCAN) IF (ISUBSCAN.GT.1) $ IFIRST_REC = IFIRST_REC + NSUBREC(ISUBSCAN-1) * * --- Raz LST_ROW,EL_ROW,WEIGHT DO I=1,NCOL LST_ROW(I) = 0 EL_ROW(I) = 0 WEIGHT_ROW(I) = 0 ENDDO * C C Add contribution of each point of subscan to the pixels of the new grid C This contribution is multiplied by a weight depending on the distance C DO I = IFIRST_REC,IFIRST_REC+NP_SUBSCAN-1 * * compute pos = coordinates of subscan point in the new grid POSX = (XSCAN_COORD(I)/10. - XMIN)/XINC + 1.D0 * IFIRST = POSX - TOTWIDX+1 ! first point of new grid to consider ILAST = POSX + TOTWIDX ! last point of new grid to consider * IF (IFIRST.LT.1) IFIRST = 1 IF (ILAST .GT.NCOL) ILAST = NCOL * DO J = IFIRST, ILAST X = FLOAT(J) - POSX ! distance between "j" and subscan poi X = ABS(X)/FUNWIDX ! convert in function unit IFUNC = INT(50. * X + 0.5) + 1 IF (IFUNC.GT.200) IFUNC = 200 Z = FUNC(IFUNC) * LST_ROW(J) = LST_ROW(J) + LST_RECORD(I)*Z EL_ROW(J) = EL_ROW(J) + EL_RECORD(I)*Z WEIGHT_ROW(J) = WEIGHT_ROW(J) + Z ! weight ENDDO ENDDO * * IF (NROW.NE.NSUBSCAN) THEN Y = (YSCAN_COORD(IFIRST_REC)/10.-YMIN)/YINC IF ( Y-NINT(Y) .LT. 0.2) THEN! test si yscan_coord du point proche POSY = NINT(Y) +1.0 ! d'une ligne de la carte (Az,EL) ELSE POSY = Y +1.0 ENDIF IFIRST = NINT(POSY) - TOTWIDY ILAST = NINT(POSY) + TOTWIDY IF (IFIRST.LT.1) IFIRST = 1 IF (ILAST.GT.NROW) ILAST = NROW DO I=1,NCOL IJ = I + (IFIRST-1)*NCOL * DO J = IFIRST, ILAST Y = (FLOAT(J) - POSY) / FUNWIDY Y = ABS(Y) IFUNC = INT(50.0 * Y + 0.5) + 1 IF (IFUNC.GT.200) IFUNC = 200 Z = FUNC(IFUNC) LST_MAP(IJ) = LST_MAP(IJ) + LST_ROW(I)*Z EL_MAP(IJ) = EL_MAP(IJ) + EL_ROW(I)*Z WEIGHT_MAP(IJ) = WEIGHT_MAP(IJ) + WEIGHT_ROW(I)*Z IJ = IJ + NCOL ENDDO ENDDO ELSE DO I=1,NCOL IJ = (ISUBSCAN-1)*NCOL LST_MAP(IJ+I) = LST_ROW(I) EL_MAP(IJ+I) = EL_ROW(I) WEIGHT_MAP(IJ+I) = WEIGHT_ROW(I) ENDDO ENDIF ENDDO ! do isubscan=1,nsubscan * * Divide by weight * DO I=1,NROW DO J=1,NCOL IJ = J+(I-1)*NCOL IF (ABS(WEIGHT_MAP(IJ)).GT.1D-10) THEN LST_MAP(IJ) = LST_MAP(IJ)/WEIGHT_MAP(IJ) EL_MAP(IJ) = EL_MAP(IJ) /WEIGHT_MAP(IJ) ELSE c write(6,*)'blanking',i,j,lst_map(ij),el_map(ij), c $ WEIGHT_MAP(ij) LST_MAP(IJ) = BLANKING EL_MAP(IJ) = BLANKING ENDIF ENDDO ENDDO RETURN END