* Last processed by NICE on 12-Jun-2000 15:51:00 * Customized for : IEEE, LINUX, UNIX, MOTIF, F77 PROGRAM MAP C------------------------------------------------------------------------ C TASK Compute a MAP from a Single-Dish Table by interpolation C and extension to 0 C C Input : C a single-dish table, C Output : C a Map C------------------------------------------------------------------------ * Global variables: INCLUDE 'inc:xpar.inc' INCLUDE 'inc:ypar.inc' INCLUDE 'inc:zpar.inc' INCLUDE 'inc:format.inc' INCLUDE 'inc:errcod.inc' INCLUDE 'inc:memory.inc' * * Local variables: * PI is defined with more digits than necessary to avoid losing * the last few bits in the decimal to binary conversion REAL*8 PI PARAMETER (PI=3.14159265358979323846D0) INTEGER SIC_GETVM, LENC EXTERNAL DUMMY * CHARACTER*80 TABLE,MAP_NAME,NAME CHARACTER*32 CHRA,CHDE REAL SD_BEAM,XMIN,XMAX,YMIN,YMAX,SEC REAL*8 NEW(2),OLD(2),XCONV(3),YCONV(3) INTEGER*4 ADDR,POINTER,IPX,IPY,IPXC,IPYC,IPW,IPR,IPI INTEGER XCOL,YCOL,WCOL,MCOL(2),LCOL,OCOL INTEGER I,IER,N,NC,NX,NY,ND,NP,NDP,MORE,NCP INTEGER*4 BLC(4), TRC(4) LOGICAL ERROR, LOGRAN(7) * DATA BLC/4*0/, TRC/4*0/ *------------------------------------------------------------------------ * Code: CALL GILDAS_OPEN CALL GILDAS_CHAR('TABLE$',TABLE) CALL GILDAS_CHAR('MAP$',MAP_NAME) CALL GILDAS_INTE('XCOL$',XCOL,1) CALL GILDAS_INTE('YCOL$',YCOL,1) CALL GILDAS_INTE('WCOL$',WCOL,1) CALL GILDAS_INTE('MCOL$',MCOL,2) CALL GILDAS_REAL('SD_BEAM$', SD_BEAM,1) ! single dish beam (rad) * CALL GILDAS_CHAR('MOSAIC_RA$',CHRA) CALL GILDAS_CHAR('MOSAIC_DEC$',CHDE) CALL GILDAS_CLOSE * * Input file N = LENC(TABLE) IF (N.LE.0) GOTO 999 NAME = TABLE(1:N) CALL SIC_PARSEF(NAME,Y_FILE,' ','.tab') CALL GDF_GEIS (Y_ISLO,ERROR) * * Cheating with Read/Write status IF (.NOT.ERROR) CALL $GDF_REIS (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR) IF (ERROR) THEN WRITE(6,*) 'F-UV_SINGLE, Cannot read input table' GOTO 999 ENDIF IF (Y_FORM.NE.FMT_R4) THEN WRITE(6,*) 'F-UV_SINGLE, Only real format supported' GOTO 999 ENDIF * * Read header CALL GDF_READY(Y_ISLO,ERROR) IF (XCOL.GT.Y_DIM(1) .OR. YCOL.GT.Y_DIM(1)) THEN WRITE(6,*) 'F-UV_SINGLE, X or Y column does not exist' GOTO 999 ENDIF ND = Y_DIM(1) NP = Y_DIM(2) CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR) * * Add a few more points on map edge with Zero values and Weights NDP = NP+8 * * Recompute offsets if needed IF (SIC_GETVM(3*NDP,ADDR).NE.1) GOTO 999 IPXC = POINTER(ADDR,MEMORY) IPYC = IPXC+NDP IPW = IPYC+NDP IPY = POINTER(Y_ADDR,MEMORY) CALL SIC_DECODE(CHRA,NEW(1),24,ERROR) IF (ERROR) GOTO 999 CALL SIC_DECODE(CHDE,NEW(2),360,ERROR) IF (ERROR) GOTO 999 OLD(1) = Y_A0 OLD(2) = Y_D0 CALL DOPOINT (MEMORY(IPY),ND,NP,XCOL,YCOL, $OLD,NEW,MEMORY(IPXC),MEMORY(IPYC)) Y_A0 = NEW(1) Y_D0 = NEW(2) * * Compute extrema and number of pixels CALL FINMAX (MEMORY(IPXC),NP,XMIN,XMAX) CALL FINMAX (MEMORY(IPYC),NP,YMIN,YMAX) SEC = PI/180.0/3600.0 XCONV(3) = -SEC YCONV(3) = +SEC * * Find size NX = 2 * MAX (NINT(ABS(XMAX/XCONV(3))+1), $NINT(ABS(XMIN/XCONV(3))+1) ) NY = 2 * MAX (NINT(ABS(YMAX/YCONV(3))+1), $NINT(ABS(YMIN/YCONV(3))+1) ) MORE = NINT(4*SD_BEAM/SEC) NX = NX+MORE NY = NY+MORE * * Extend to nearest power of two I = 32 DO WHILE(I.LT.NX) I = I*2 ENDDO NX = I I = 32 DO WHILE(I.LT.NY) I = I*2 ENDDO NY = I XCONV(1) = NX/2+1 XCONV(2) = 0.0 YCONV(1) = NY/2+1 YCONV(2) = 0.0 * * Define coordinates and values of additional data points CALL SETMORE (MEMORY(IPXC),MEMORY(IPYC),NP,NDP, $NX,XCONV,NY,YCONV) WRITE(6,*) 'Creating a cube with ',NX,' by ',NY,' pixels' WRITE(6,*) 'Increment ',0.1*NINT(YCONV(3)*10*180*3600/PI) * * Warn for big images IF (NX.GT.4096 .OR. NY.GE.4096) THEN WRITE(6,*) 'E-UV_SINGLE, More than 4096 pixels in X or Y' GOTO 999 ELSEIF (NX.GT.512 .OR. NY.GT.512) THEN WRITE(6,*) 'W-UV_SINGLE, More than 512 pixels in X or Y' ENDIF * * Define output channels. One more for weights IF (MCOL(2).EQ.0) MCOL(2) = ND MCOL(1) = MAX(1,MIN(MCOL(1),ND)) MCOL(2) = MAX(1,MIN(MCOL(2),ND)) OCOL = MIN(MCOL(1),MCOL(2))-1 LCOL = MAX(MCOL(1),MCOL(2)) CALL GDF_CHYX NC = LCOL-OCOL+1 WRITE(6,*) 'Creating ',NC,' Channels from ',MCOL * * Create image (in order L M V ) NAME = MAP_NAME(1:N) CALL SIC_PARSEF(NAME,X_FILE,' ','.lmv') WRITE(6,*) 'I-UV_SINGLE, Creating map file '// $X_FILE(1:LENC(X_FILE)) X_NDIM = 3 X_DIM(1) = NX X_DIM(2) = NY X_DIM(3) = NC X_DIM(4) = 1 X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3) X_REF3 = Y_REF1-OCOL X_VAL3 = X_VOFF X_INC3 = X_VRES X_REF1 = XCONV(1) X_VAL1 = XCONV(2) X_INC1 = XCONV(3) X_REF2 = YCONV(1) X_VAL2 = YCONV(2) X_INC2 = YCONV(3) X_CODE(1) = Y_CODE(2) X_CODE(2) = Y_CODE(3) X_CODE(3) = Y_CODE(1) X_GENE = 29 ! Not a table stupid... X_EXTR = 0 ! extrema not computed X_XAXI = 1 ! Reset projected axis X_YAXI = 2 ! X_FAXI = 3 * CALL GDF_GEIS (X_ISLO,ERROR) CALL GDF_WRITX(X_ISLO,ERROR) X_FORM = FMT_R4 X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4) IF (.NOT.ERROR) CALL $GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR) IF (ERROR) THEN WRITE(6,*) 'F-UV_SINGLE, Cannot create output LMV image' GOTO 999 ENDIF * * Compute weights plane and triangulation BLC(3) = NC BLC(4) = 1 TRC(3) = NC TRC(4) = 1 CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR) IPX = POINTER(X_ADDR,MEMORY) * * Load weights CALL DODATA (MEMORY(IPY),ND,NP,WCOL,MEMORY(IPW),NDP) * * Initialize gridding IER = SIC_GETVM (13*NDP,ADDR) IF (IER.NE.1) GOTO 999 IPR = POINTER(ADDR,MEMORY) IER = SIC_GETVM (35*NDP,ADDR) IF (IER.NE.1) GOTO 999 IPI = POINTER(ADDR,MEMORY) CALL GRIDINI(NX,XCONV(1),XCONV(2),XCONV(3), $NY,YCONV(1),YCONV(2),YCONV(3)) LOGRAN(7) = .TRUE. ! Extrapolate if needed LOGRAN(3) = .FALSE. ! Compute triangulation NCP = 4 CALL GRIDSET(NCP,X_BVAL,LOGRAN) CALL GRIDRAN(MEMORY(IPXC), $MEMORY(IPYC), $MEMORY(IPW), $NDP, $MEMORY(IPR), $MEMORY(IPI), $MEMORY(IPX), $DUMMY,DUMMY,DUMMY,ERROR) LOGRAN(3) = .TRUE. ! Just interpolate for others CALL GDF_FRMS(X_MSLO,ERROR) * DO I=1,NC-1 BLC(3) = I TRC(3) = I CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR) IPX = POINTER(X_ADDR,MEMORY) * * Load data values CALL DODATA (MEMORY(IPY),ND,NP,OCOL+I-1,MEMORY(IPW),NDP) CALL GRIDRAN(MEMORY(IPXC), $ MEMORY(IPYC), $ MEMORY(IPW), $ NDP, $ MEMORY(IPR), $ MEMORY(IPI), $ MEMORY(IPX), $ DUMMY,DUMMY,DUMMY,ERROR) CALL GDF_FRMS(X_MSLO,ERROR) ENDDO CALL GDF_FRIS(X_ISLO,ERROR) CALL GDF_FRIS(Y_ISLO,ERROR) CALL GAGOUT('I-MAP, Successful completion') * 999 CALL GAGOUT('E-MAP, Fatal error') CALL SYSEXI(FATALE) END * SUBROUTINE DOPOINT(INPUT,ND,NP,XCOL,YCOL,OLD,NEW,XD,YD) INTEGER ND,NP,XCOL,YCOL REAL INPUT(ND,NP),XD(NP),YD(NP) REAL*8 OLD(2),NEW(2) * REAL*8 DRA,DDE,RA,DE,UNCDE,CDE INTEGER I * UNCDE = 1.D0/COS(OLD(2)) CDE = COS(NEW(2)) * DO I=1,NP RA = OLD(1) + DBLE(INPUT(XCOL,I))*UNCDE DE = OLD(2) + DBLE(INPUT(YCOL,I)) DRA = (RA - NEW(1)) * CDE DDE = DE - NEW(2) XD(I) = DRA YD(I) = DDE ENDDO END * SUBROUTINE SETMORE(XD,YD,NP,NDP,NX,XCONV,NY,YCONV) INTEGER NP,NDP,NX,NY REAL XD(NDP),YD(NDP) REAL*8 XCONV(3),YCONV(3) * XD(NP+1) = (1.0-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+2) = (NX/2-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+3) = (NX-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+4) = (1.0-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+5) = (NX-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+6) = (1.0-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+7) = (NX/2-XCONV(1))*XCONV(3)+XCONV(2) XD(NP+8) = (NX-XCONV(1))*XCONV(3)+XCONV(2) * YD(NP+1) = (1.0-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+2) = (1.0-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+3) = (1.0-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+4) = (NY/2-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+5) = (NY/2-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+6) = (NY-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+7) = (NY-YCONV(1))*YCONV(3)+YCONV(2) YD(NP+8) = (NY-YCONV(1))*YCONV(3)+YCONV(2) * END * SUBROUTINE DODATA (INPUT,ND,NP,IC,OUTPUT,NDP) INTEGER ND,NP,NDP,IC REAL INPUT(ND,NP) REAL OUTPUT (NDP) * INTEGER I * DO I=1,NP OUTPUT(I) = INPUT(IC,I) ENDDO DO I=NP+1,NDP OUTPUT(I) = 0.0 ENDDO END * SUBROUTINE DUMMY END * SUBROUTINE FINMAX (X,NP,XMIN,XMAX) INTEGER NP,I REAL X(NP),XMIN,XMAX * XMIN = X(1) XMAX = X(1) DO I=2,NP IF (X(I).LT.XMIN) THEN XMIN = X(I) ELSEIF (X(I).GT.XMAX) THEN XMAX = X(I) ENDIF ENDDO END