LOCAL INCLUDE 'EVAUV.INC' C Include EVAUV INCLUDE 'INCS:PUVD.INC' INCLUDE 'INCS:ZPBUFSZ.INC' INTEGER UVHIST, MAHIST PARAMETER (UVHIST = 302) PARAMETER (MAHIST = 402) C NOTE: uses PARAMETER in DGDS.INC REAL XSIN, XDISIN, XQUAL, XTIME(8), XBAND, XFREQ, XFQID, * XANT(50), XBASE(50), XSUBA, XBIF, XEIF, XBCHAN, XECHAN, * XDOCAL, XGUSE, XDOPOL, XPDVER, XBLVER, XFLAG, XDOBND, XBPVER, * XSMOTH(3), XNMAPS, XNGAUS, XS2, XDISK2, XVER, XBCOMP(MAXAFL), * XNCOMP(MAXAFL), XFLUX, SMODEL(7), SPIX(4), DOOUT, XSO, XDISKO, * DOHIST, XCELL(2), APARM(10), PXRANG(2), XLABEL, XDOTV, * XGRCH, BADD(10), XTRPRM(10), RPARM(30) HOLLERITH XNAMEI(3), XCLAIN(2), XSOUR(4,30), XCALC, XXSTOK, * XNAME2(3), XCLAS2(2), XCMETH, XOUTNM(3), XFUNC CHARACTER NAMEIN*12, CLAIN*6, NAME2*12, CLAS2*6, CMETH*4, CMOD*4, * OUTNAM*12, FUNCTY*2, MUNITS*8 REAL BUFF1(UVBFSS), BUFF2(UVBFSS), BUFF3(UVBFSS), * BUFF4(UVBFSS), HISTA(UVHIST,2), HISTR(UVHIST,2), RSCALE, * MSCALE(2,2), MFIT(2), REMXN(2,2), IMMXN(2,2) INTEGER NCOMP(MAXFLD), BCOMP(MAXFLD), UVBLK(256), DISKIN, DISK2, * OLDCNO, CLNCNO, JBUFSZ, SEQIN, SEQ2, VER, CHAN, NCHAN, ISTOKE, * METHOD, MODEL, LBIF, LEIF, CH1, ISLFRQ, NGAUSS, SC1FIL, SC2FIL, * HISTM(MAHIST,2), NPARM, IBUFF1(UVBFSS), IBUFF4(UVBFSS), WVOL, * WCNO, WSCI, INCATB(256) DOUBLE PRECISION INCATD(128) EQUIVALENCE (INCATD, INCATB) EQUIVALENCE (IBUFF1, BUFF1), (IBUFF4, BUFF4) COMMON /INFO/ INCATB, UVBLK, BCOMP, NCOMP, DISKIN, DISK2, OLDCNO, * CLNCNO, ISTOKE, JBUFSZ, SEQIN, SEQ2, VER, CHAN, NCHAN, METHOD, * MODEL, LBIF, LEIF, CH1, ISLFRQ, NGAUSS, SC1FIL, SC2FIL, HISTA, * HISTR, HISTM, RSCALE, MSCALE, NPARM, MFIT, REMXN, IMMXN, WVOL, * WCNO, WSCI COMMON /BUFRS/ BUFF1, BUFF2, BUFF3, BUFF4 COMMON /INPARM/ XNAMEI, XCLAIN, XSIN, XDISIN, XSOUR, XQUAL, XCALC, * XXSTOK, XTIME, XBAND, XFREQ, XFQID, XANT, XBASE, XSUBA, XBIF, * XEIF, XBCHAN, XECHAN, XDOCAL, XGUSE, XDOPOL, XPDVER, XBLVER, * XFLAG, XDOBND, XBPVER, XSMOTH, XNMAPS, XNGAUS, * XNAME2, XCLAS2, XS2, XDISK2, XVER, XBCOMP, XNCOMP, XFLUX, * XCMETH, SMODEL, SPIX, DOOUT, XOUTNM, XSO, XDISKO, DOHIST, * XCELL, APARM, PXRANG, XFUNC, XLABEL, XDOTV, XGRCH, BADD, * XTRPRM, RPARM COMMON /CHRCOM/ NAMEIN, CLAIN, NAME2, CLAS2, CMETH, CMOD, MUNITS, * OUTNAM, FUNCTY C End EVAUV LOCAL END LOCAL INCLUDE 'DSPX.INC' C include for spectral index C used with SMODEL REAL PTSPIX(4) LOGICAL DOPTSP COMMON /PTSPEC/ PTSPIX, DOPTSP LOCAL END PROGRAM EVAUV C----------------------------------------------------------------------- C! Subtracts & divides a model into uv data - evaluates statistics C# UV AP-appl Modeling C----------------------------------------------------------------------- C; Copyright (C) 2009-2012, 2014-2015, 2017-2020 C; Associated Universities, Inc. Washington DC, USA. C; C; This program is free software; you can redistribute it and/or C; modify it under the terms of the GNU General Public License as C; published by the Free Software Foundation; either version 2 of C; the License, or (at your option) any later version. C; C; This program is distributed in the hope that it will be useful, C; but WITHOUT ANY WARRANTY; without even the implied warranty of C; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C; GNU General Public License for more details. C; C; You should have received a copy of the GNU General Public C; License along with this program; if not, write to the Free C; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C; MA 02139, USA. C; C; Correspondence concerning AIPS should be addressed as follows: C; Internet email: aipsmail@nrao.edu. C; Postal address: AIPS Project Office C; National Radio Astronomy Observatory C; 520 Edgemont Road C; Charlottesville, VA 22903-2475 USA C----------------------------------------------------------------------- C EVAUV is a task to subtract a source model from and divide a source C model into a uv-data base. It then does statistics on the results. C It also does statistics on the model images. C Inputs: C AIPS adverb Prg. name. Description. C INNAME NAMEIN Name of input UV data. C INCLASS CLASIN Class of input UV data. C INSEQ SEQIN Seq. of input UV data. C INDISK DISKIN Disk number of input UV data. C NMAPS MFIELD Number of input images. C NGAUSS NGAUS Number Gaussians in NMAPS C IN2NAME NAME2 Name of map with CLEAN components. C IN2CLASS CLAS2 Class of map with CLEAN components. C IN2SEQ SEQ2 Seq. of map with CLEAN components. C IN2DISK DISK2 Vol. of map with CLEAN components. C INVER VER Version no. of CC file. C BCOMP(64) BCOMP Start clean component to sub. C 1 per field. C NCOMP(64) NCOMP Last Clean component no to sub. C 1 per field, 0 => all C CMETHOD METHOD Modeling method: C 'DFT' = FDT method C 'GRID' = gridded FFT method. C ' ' chose fastest. C CMODEL MODEL Model type, 'COMP'=>CC C 'IMAG'=> image. C BADDISK IBAD Disk nos. to avoid for scratch files C----------------------------------------------------------------------- CHARACTER PRGM*6, UTYPE*2 INTEGER IRET, IERR, N, DSK1, DSK2, CNO1, CNO2 LOGICAL DOMSG DOUBLE PRECISION DMAP(1025,1025), APCORE(2) REAL SCRIMG(1025,1025,2) INCLUDE 'INCS:PUVD.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'EVAUV.INC' INCLUDE 'DSPX.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DSCD.INC' INCLUDE 'INCS:DTVC.INC' INCLUDE 'INCS:DCAT.INC' REAL BUFF5(MAXIF) DATA PRGM /'EVAUV '/ DATA DOMSG /.TRUE./ C----------------------------------------------------------------------- C Get input parameters and C create output scratch files CALL EVAUVI (PRGM, IRET) IF (IRET.NE.0) THEN RQUICK = .FALSE. GO TO 990 END IF DATDIV = .TRUE. IF (SC1FIL.GT.0) THEN DSK1 = 0 CNO1 = SC1FIL DSK2 = 0 CNO2 = SC2FIL ELSE DSK1 = FVOL(-SC1FIL) CNO1 = FCNO(-SC1FIL) DSK2 = FVOL(-SC2FIL) CNO2 = FCNO(-SC2FIL) END IF C Process data CALL UVMDVS (APCORE, -WVOL, WSCI, DSK1, CNO1, DSK2, CNO2, MODEL, * METHOD, DOMSG, CHAN, NCHAN, JBUFSZ, ISLFRQ, BUFF1, BUFF2, * BUFF3, IBUFF4, BUFF5, IRET) C do image stats: get UV size IF ((IRET.EQ.0) .AND. (.NOT.DOPTMD)) CALL EVAIMG (NGAUSS, MAHIST, * PXRANG, MSCALE, MUNITS, HISTM, MFIT, RPARM, BUFF1, IRET) IF (.NOT.DOPTMD) CALL UNSETG (BUFF2) C do UV stats IF (IRET.EQ.0) CALL EVAUVS (IRET) C do plots IF (IRET.EQ.0) THEN IF (RQUICK) THEN CALL PTPARM (30, RPARM, BUFF1, IERR) CALL RELPOP (IRET, BUFF1, IERR) END IF N = APARM(3) + 0.5 N = MIN (N, 1025) N = (N / 2) * 2 + 1 CALL EVAPLT (N, DMAP, SCRIMG(1,1,1), SCRIMG(1,1,2), IRET) IRET = MAX (0, IRET) C Put input file in READ IF (DOHIST.GT.0.0) THEN FRW(1) = 0 UTYPE = 'UV' CALL CATDIR ('CSTA', DISKIN, OLDCNO, NAMEIN, CLAIN, SEQIN, * UTYPE, NLUSER, 'CLWR', BUFF1, IERR) CALL CATDIR ('CSTA', DISKIN, OLDCNO, NAMEIN, CLAIN, SEQIN, * UTYPE, NLUSER, 'READ', BUFF1, IERR) END IF ELSE RQUICK = .FALSE. END IF C HI files etc IF ((IRET.EQ.0) .AND. (DOOUT.GT.0.0)) CALL EVAUVH C Close down 990 IF (.NOT.RQUICK) CALL PTPARM (30, RPARM, BUFF1, IERR) CALL DIE (IRET, BUFF1) C 999 STOP END SUBROUTINE EVAUVI (PRGN, IERR) C----------------------------------------------------------------------- C EVAUVI gets input parameters for EVAUV and creates output scratch C files. C Inputs: C PRGN C*6 Task name C Output: C IERR I Error code: non-zero => quit C----------------------------------------------------------------------- CHARACTER PRGN*6 INTEGER IERR C CHARACTER STAT*4, UTYPE*2 INTEGER JERR, INMETH, I, IBUFF(512), IROUND, MXFLD, ISIZE, * DISKO, SEQOUT LOGICAL T, F, WASOME DOUBLE PRECISION DEMP INCLUDE 'INCS:PUVD.INC' INCLUDE 'EVAUV.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DCAT.INC' INCLUDE 'DSPX.INC' EQUIVALENCE (IBUFF, BUFF2) DATA T, F /.TRUE.,.FALSE./ C----------------------------------------------------------------------- C Initialize I/O CALL ZDCHIN (T, BUFF1) CALL VHDRIN JBUFSZ = UVBFSS * 2 MXFLD = MAXAFL ISLFRQ = -1 C Initialize /CFILES/ NSCR = 0 NCFILE = 0 IERR = 0 JERR = 0 C Get input parameters. NPARM = 61 + 2 * MXFLD + 254 CALL GTPARM (PRGN, NPARM, RQUICK, XNAMEI, BUFF1, IERR) IF (IERR.NE.0) THEN RQUICK = .FALSE. JERR = 8 IF (IERR.EQ.1) GO TO 999 WRITE (MSGTXT,1000) IERR CALL MSGWRT (8) END IF C Restart AIPS IF (JERR.NE.0) GO TO 999 IERR = 5 CALL RFILL (30, 0.0, RPARM) C Crunch input parameters. C Convert characters CALL H2CHR (12, 1, XNAMEI, NAMEIN) CALL H2CHR (6, 1, XCLAIN, CLAIN) CALL H2CHR (12, 1, XNAME2, NAME2) CALL H2CHR (6, 1, XCLAS2, CLAS2) CALL H2CHR (4, 1, XCMETH, CMETH) CMOD = ' ' CALL H2CHR (2, 1, XFUNC, FUNCTY) CALL H2CHR (12, 1, XOUTNM, OUTNAM) C File sequence numbers SEQIN = IROUND (XSIN) SEQ2 = IROUND (XS2) C File disk numbers DISKIN = IROUND (XDISIN) DISK2 = IROUND (XDISK2) C CC file version number VER = IROUND (XVER) C Default LEIF when have CATBLK C Number of fields MFIELD = 1 IF (XNMAPS.GT.0.0) MFIELD = IROUND (XNMAPS) NGAUSS = XNGAUS NGAUSS = MAX (1, NGAUSS) LIMFLX = XFLUX C Start component number NONEG = F WASOME = F DO 8 I = 1,MFIELD IF (I.LE.MAXAFL) THEN BCOMP(I) = XBCOMP(I) + 0.1 BCOMP(I) = MAX (BCOMP(I), 1) NCOMP(I) = ABS (XNCOMP(I)) + 0.1 IF (XNCOMP(I).LE.-0.5) NONEG = T IF (NCOMP(I).GT.0) WASOME = T ELSE BCOMP(I) = 1 NCOMP(I) = 0 IF (WASOME) NCOMP(I) = 1000000000 END IF 8 CONTINUE C Factor FACGRD(1) = 1.0 FACGRD(2) = 1.0 C Disks to avoid for scratch DO 10 I = 1,10 IBAD(I) = IROUND (BADD(I)) 10 CONTINUE C Get modeling method METHOD = 0 IF (CMETH.EQ.'DFT ') METHOD = -1 IF (CMETH.EQ.'GRID') METHOD = 1 C Get model type MODEL = 0 IF (CMOD.EQ.'COMP') MODEL = 1 IF (CMOD.EQ.'IMAG') MODEL = 2 IF (CMOD.EQ.'SUBI') MODEL = 3 C Get CATBLK from old file. OLDCNO = 1 UTYPE = 'UV' CALL CATDIR ('SRCH', DISKIN, OLDCNO, NAMEIN, CLAIN, SEQIN, UTYPE, * NLUSER, STAT, BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1010) IERR, NAMEIN, CLAIN, SEQIN, 'UV', * DISKIN, NLUSER GO TO 990 END IF CALL CATIO ('READ', DISKIN, OLDCNO, CATBLK, 'REST', BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1020) IERR GO TO 990 END IF NCFILE = NCFILE + 1 FVOL(NCFILE) = DISKIN FCNO(NCFILE) = OLDCNO FRW(NCFILE) = 0 CALL COPY (256, CATBLK, INCATB) C Get uv header info. CALL UVPGET (IERR) IF (IERR.NE.0) GO TO 999 C copy file with cal/flag/etc CALL EVCOPY (IERR) IF (IERR.NE.0) GO TO 999 C Reject multisource files (it is C not strictly necessary to do C this but input files can not C contain more than one source C and EVAUV would need to verify C this here) IF (ILOCSU.GE.0) THEN MSGTXT = 'INPUT FILE MUST BE A SINGLE-SOURCE FILE.' CALL MSGWRT (8) MSGTXT = 'USE SPLIT TO EXTRACT DESIRED DATA.' CALL MSGWRT (8) IERR = 1 GO TO 999 END IF C C channel number CHAN = 1 NCHAN = CATBLK(KINAX+JLOCF) C IFs LBIF = 1 LEIF = 1 C Deal with IFs IF (JLOCIF.GE.0) THEN LEIF = CATBLK(KINAX+JLOCIF) C Offset CHAN for IF IF (JLOCIF.LT.JLOCF) THEN C do all IFs if all ch in this C axis order IF (NCHAN.GT.1) THEN LBIF = 1 LEIF = CATBLK(KINAX+JLOCIF) END IF CH1 = CATBLK(KINAX+JLOCIF) * (CHAN-1) + LBIF C chan first usually ELSE CH1 = CATBLK(KINAX+JLOCF) * (LBIF-1) + CHAN END IF C Change NCHAN to include IFs NCHAN = NCHAN * (LEIF - LBIF + 1) C Change Start CHAN for IFs CHAN = CH1 C Reset INCF INCF = MIN (INCF, INCIF) END IF C Check order of u,v,w IF (((ILOCV-ILOCU).NE.1) .OR. ((ILOCW-ILOCV).NE.1)) THEN WRITE (MSGTXT,1070) ILOCU, ILOCV, ILOCW IERR = 1 GO TO 990 END IF C Setup common for modeling C routines C Uv header block CALL COPY (256, CATBLK, UVBLK) INMETH = METHOD IF (APARM(1).LE.0.0) APARM(1) = 5.0 IF (APARM(2).LE.0.0) APARM(2) = 10.0 IF (APARM(3).LE.0.0) THEN DEMP = NCHAN DEMP = DEMP * NVIS DEMP = SQRT (DEMP/10.0D0) DEMP = MAX (50.0D0, MIN (1024.0D0, DEMP)) APARM(3) = DEMP END IF I = APARM(4) + 0.1 I = MIN (I, 33) I = (I / 2) * 2 + 1 APARM(4) = I C Make sure method='GRID' for C image model. IF (MODEL.GE.2) THEN MSGTXT = 'Warning: GRID used for Images models' IF (METHOD.EQ.-1) CALL MSGWRT(5) METHOD = 1 END IF C Check model DOPTMD = ABS (SMODEL(1)) .GT. 1.0E-20 PTFLX = SMODEL(1) PTRAOF = SMODEL(2) PTDCOF = SMODEL(3) PARMOD(1) = SMODEL(4) PARMOD(2) = SMODEL(5) PARMOD(3) = SMODEL(6) PARMOD(4) = SMODEL(7) IF (DOPTMD) THEN CALL RCOPY (4, SPIX, PTSPIX) C get spectrum IF (PTFLX.LT.0.0) THEN CALL GETSPX (PTFLX, PTSPIX, IERR) IF (IERR.NE.0) GO TO 999 SMODEL(1) = PTFLX CALL COPY (2, PTSPIX, SPIX) END IF DOPTSP = (SPIX(1).NE.0.0) .OR. (SPIX(2).NE.0.0) RPARM(26) = SMODEL(1) CALL RCOPY (4, SPIX, RPARM(27)) ELSE CALL RFILL (4, 0.0, PTSPIX) CALL RFILL (5, 0.0, RPARM(26)) DOPTSP = .FALSE. END IF ISTOKE = 1 C If NONEG or DOPTMD use DFT IF (DOPTMD) THEN METHOD = -1 WRITE (MSGTXT,1090) SMODEL(1), SMODEL(2), SMODEL(3) CALL MSGWRT (4) MSGTXT = 'IGNORING MODEL IMAGE ' // NAME2 // CLAS2 IF ((NAME2.NE.' ') .OR. (CLAS2.NE.' ')) CALL MSGWRT (7) IF (DOPTSP) THEN WRITE (MSGTXT,1091) PTSPIX CALL MSGWRT (4) END IF END IF C Warn user if changing METHOD MSGTXT = '*** WARNING: OVERRIDING SPECIFIED COMPUTATION METHOD' IF ((INMETH.NE.METHOD) .AND. (INMETH.NE.0)) CALL MSGWRT (6) C Get info on model file(s) IF (DOPTMD) THEN DO3DIM = .FALSE. ELSE CALL SETGDS (WVOL, WCNO, NAME2, CLAS2, SEQ2, DISK2, MFIELD, * VER, NCOMP, BCOMP, MODEL, METHOD, BUFF1, BUFF2, ISTOKE, * IERR) IF (IERR.NE.0) GO TO 999 IF (MODEL.EQ.3) THEN MSGTXT = 'Using sub-images for the source model' ELSE IF (MODEL.EQ.2) THEN MSGTXT = 'Using images for the source model' ELSE MSGTXT = 'Using Clean Component source model' END IF CALL MSGWRT (3) CALL FACSET (WVOL, WCNO, 1, 0, MODEL, FACGRD, IERR) IF (IERR.NE.0) GO TO 999 END IF C Create scratch files IF (DOOUT.LE.0.0) THEN CALL UVSIZE (LREC, NVIS, ISIZE) CALL SCREAT (ISIZE, BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1130) IERR, 'SCRATCH' GO TO 990 END IF SC1FIL = NSCR CALL SCREAT (ISIZE, BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1130) IERR, 'SCRATCH' GO TO 990 END IF SC2FIL = NSCR C cataloged files out ELSE CALL COPY (256, UVBLK, CATBLK) SEQOUT = IROUND (XSO) DISKO = IROUND (XDISKO) IF (OUTNAM.EQ.' ') OUTNAM = NAMEIN CALL CHR2H (12, OUTNAM, KHIMNO, CATH(KHIMN)) CALL CHR2H (6, 'EVAUVS', KHIMCO, CATH(KHIMC)) CATBLK(KIIMS) = SEQOUT CCNO = 1 CALL UVCREA (DISKO, CCNO, BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1130) IERR, 'CATALOGED' GO TO 990 END IF NCFILE = NCFILE + 1 FVOL(NCFILE) = DISKO FCNO(NCFILE) = CCNO FRW(NCFILE) = 2 SC1FIL = -NCFILE C copy keywords CALL KEYCOP (WVOL, WCNO, DISKO, CCNO, IERR) C division file CALL COPY (256, UVBLK, CATBLK) SEQOUT = IROUND (XSO) DISKO = IROUND (XDISKO) CALL CHR2H (12, OUTNAM, KHIMNO, CATH(KHIMN)) CALL CHR2H (6, 'EVAUVD', KHIMCO, CATH(KHIMC)) CATBLK(KIIMS) = SEQOUT CCNO = 1 CALL UVCREA (DISKO, CCNO, BUFF1, IERR) IF (IERR.NE.0) THEN WRITE (MSGTXT,1130) IERR, 'CATALOGED' GO TO 990 END IF NCFILE = NCFILE + 1 FVOL(NCFILE) = DISKO FCNO(NCFILE) = CCNO FRW(NCFILE) = 2 SC2FIL = -NCFILE C copy keywords CALL KEYCOP (WVOL, WCNO, DISKO, CCNO, IERR) CALL COPY (256, UVBLK, CATBLK) END IF C Put input file in READ IF (DOHIST.GT.0.0) THEN STAT = 'WRIT' FRW(1) = 1 ELSE STAT = 'READ' END IF UTYPE = 'UV' CALL CATDIR ('CSTA', DISKIN, OLDCNO, NAMEIN, CLAIN, SEQIN, UTYPE, * NLUSER, STAT, BUFF1, IERR) IERR = 0 GO TO 999 C Error messages output 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVAUVI: ERROR',I3,' OBTAINING INPUT PARAMETERS') 1010 FORMAT ('ERROR',I3,' FINDING ',A12,' . ',A6,' . ', * I3,2X,A2,' DISK=',I3,' USID=',I4) 1020 FORMAT ('ERROR',I3,' COPYING CATBLK ') 1070 FORMAT ('WRONG ORDER FOR U, V, W =',3I4) 1090 FORMAT ('Using SMODEL =',F10.5,2F12.5) 1091 FORMAT ('With spectral index =',4F9.4) 1130 FORMAT ('ERROR',I3,' CREATING OUTPUT ',A,' FILE') END SUBROUTINE GETSPX (PTFLX, PTSPIX, IERR) C----------------------------------------------------------------------- C GETSPX tries to get the flux and spectral index parameters from the C SU table C In/out: C PTFLX R Flux at obs frequency: -2 in -> curvature C Output C PTSPIX R(2) Spectral index, curvature C IERR I Failed C----------------------------------------------------------------------- INTEGER IERR REAL PTFLX, PTSPIX(*) C INCLUDE 'INCS:DSEL.INC' REAL SPIX(4) INTEGER NTERM DOUBLE PRECISION DT, X INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DMSG.INC' C----------------------------------------------------------------------- NTERM = 1 IF (PTFLX.LT.-1.5) NTERM = 2 CALL FNDSPX (IUDISK, IUCNO, SOUWAN, FRQSEL, CATUV, NTERM, SPIX, * IERR) IF (FREQ.LE.0.0D0) THEN IERR = 10 MSGTXT = 'GETSPX: FREQ IS ZERO QUITTING' CALL MSGWRT (8) GO TO 999 END IF DT = LOG10 (FREQ/1.D9) X = SPIX(2) * DT + SPIX(3) * DT * DT PTFLX = SPIX(1) * (10.0D0 ** X) IF (NTERM.LE.1) THEN PTSPIX(1) = SPIX(2) PTSPIX(2) = 0.0 ELSE PTSPIX(1) = SPIX(2) + 2.0D0 * DT * SPIX(3) PTSPIX(2) = SPIX(3) END IF C 999 RETURN END SUBROUTINE EVCOPY (IRET) C----------------------------------------------------------------------- C EVCOPY interprets the input adverbs that apply data selection, C flagging, and calibration. It creates a scratch file to receive C the data sp processed and fills it. C Output: C IRET I Error code C----------------------------------------------------------------------- INTEGER IRET C INCLUDE 'EVAUV.INC' INTEGER NXANT, NXBAS, IXANT(50), IXBAS(50), LUN, NFREQ, FIND, * BO, LENBU, BIND, VO, IPTRO, NIOLIM, NIOUT, IA1, IA2, CNOSCR, * NCOPY, IROUND, ISIZE, IERR, I LOGICAL DESEL, T, F, REQBAS, MATCH REAL VPARM(20), BASEN CHARACTER NAME*48 DOUBLE PRECISION XCOUNT INCLUDE 'INCS:DSEL.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DHDR.INC' DATA VO, BO /0, 1/ DATA T, F /.TRUE.,.FALSE./ C----------------------------------------------------------------------- C Info for UVGET: CALL SELINI C Put selection criteria into C correct common. CALL H2CHR (4, 1, XCALC, SELCOD) CALL H2CHR (4, 1, XXSTOK, STOKES) DO 10 I = 1,30 CALL H2CHR (16, 1, XSOUR(1,I), SOURCS(I)) 10 CONTINUE SELQUA = IROUND (XQUAL) UNAME = NAMEIN UCLAS = CLAIN UDISK = DISKIN USEQ = SEQIN DOCAL = XDOCAL.GT.0.0 DOWTCL = DOCAL .AND. (XDOCAL.LE.99.0) C Set time range. CALL RCOPY (8, XTIME, TIMRNG) IF ((TIMRNG(1)+TIMRNG(2)+TIMRNG(3)+TIMRNG(4)) .EQ.0.0) * TIMRNG(1)=-1.0E6 IF ((TIMRNG(5)+TIMRNG(6)+TIMRNG(7)+TIMRNG(8)) .EQ.0.0) * TIMRNG(5)=1.0E6 TSTART = TIMRNG(1) + TIMRNG(2) / 24. + TIMRNG(3) / (24. * 60.) + * TIMRNG(4) / (24. * 60. * 60.) TEND = TIMRNG(5) + TIMRNG(6) / 24. + TIMRNG(7) / (24. * 60.) + * TIMRNG(8) / (24. * 60. * 60.) DOPOL = IROUND(XDOPOL) IF (XDOPOL.GT.0.0) DOPOL = MAX (1, DOPOL) PDVER = IROUND (XPDVER) DOAPPL = .FALSE. SUBARR = IROUND (XSUBA) IF (SUBARR.LT.0) SUBARR = 0 FGVER = IROUND (XFLAG) DOBAND = IROUND (XDOBND) BPVER = IROUND (XBPVER) CALL RCOPY (3, XSMOTH, SMOOTH) CLVER = IROUND (XGUSE) CLUSE = IROUND (XGUSE) BLVER = IROUND (XBLVER) C Channel selection? IF (JLOCIF.LT.0) THEN BIF = 1 EIF = 1 ELSE BIF = IROUND (XBIF) EIF = IROUND (XEIF) BIF = MIN (MAX (1, BIF), CATBLK(KINAX+JLOCIF)) IF (EIF.LT.BIF) EIF = CATBLK(KINAX+JLOCIF) END IF NFREQ = CATBLK(KINAX+JLOCF) BCHAN = IROUND (XBCHAN) ECHAN = IROUND (XECHAN) IF ((BCHAN.LE.0) .OR. (BCHAN.GT.NFREQ)) BCHAN = 1 IF ((ECHAN.LE.0) .OR. (ECHAN.GT.NFREQ)) ECHAN = NFREQ IF (BCHAN.GT.ECHAN) THEN MSGTXT = 'INVALID BCHAN AND ECHAN' CALL MSGWRT (6) IRET = 1 GO TO 990 END IF C Freq id IF (XBAND.GT.0.0) SELBAN = XBAND IF (XFREQ.GT.0.0) SELFRQ = XFREQ FRQSEL = IROUND (XFQID) IF (FRQSEL.EQ.0) FRQSEL = -1 LUN = 78 CALL FQMATC (DISKIN, OLDCNO, CATBLK, LUN, SELBAN, SELFRQ, MATCH, * FRQSEL, IRET) IF (.NOT.MATCH) THEN MSGTXT = 'NO MATCH TO SELBAND/SELFREQ ADVERBS - CHECK INPUTS' IRET = 1 GO TO 990 END IF IF (IRET.GT.0) GO TO 999 ISLFRQ = FRQSEL C Find baselines to copy CALL SETANT (50, XANT, XBASE, NXANT, NXBAS, IXANT, IXBAS, DESEL) C now using cal system - C UVGET makes header CALL UVGET ('INIT', VPARM, BUFF1, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'OPEN WITH UVGET' GO TO 990 END IF C Determine size. LENBU = 1 CALL UVSIZE (LREC, NVIS, ISIZE) C Create output file if necessary C Create scratch file. CALL SCREAT (ISIZE, BUFF1, IRET) CNOSCR = NSCR IF (IRET.NE.0) THEN IF (IRET.EQ.1) THEN MSGTXT = 'EVCOPY: TOO LITTLE DISK SPACE FOR SCRATCH FILE' ELSE WRITE (MSGTXT,1000) IRET, 'CREATING SCRATCH FILE' END IF GO TO 990 END IF C Update CATBLK: ignore error CALL CATIO ('UPDT', SCRVOL(CNOSCR), SCRCNO(CNOSCR), CATBLK, * 'REST', BUFF1, IRET) WVOL = SCRVOL(CNOSCR) WCNO = SCRCNO(CNOSCR) WSCI = CNOSCR CALL ZPHFIL ('SC', WVOL, WCNO, 1, NAME, IRET) C Open output file. CALL ZOPEN (LUN, FIND, WVOL, NAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'OPENING OUTPUT FILE' GO TO 990 END IF C Init vis file for write LENBU = 0 CALL UVINIT ('WRIT', LUN, FIND, NVIS, VO, LREC, LENBU, JBUFSZ, * BUFF2, BO, BIND, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'INIT OUTPUT FILE' GO TO 990 END IF IPTRO = BIND NIOUT = 0 NIOLIM = LENBU XCOUNT = 0.0D0 NCOPY = LREC - NRPARM C Loop C Read vis. record. 100 CALL UVGET ('READ', VPARM, BUFF1, IRET) IF (IRET.GT.0) THEN WRITE (MSGTXT,1000) IRET, 'READING VIS DATA' GO TO 990 C Loop over buffer ELSE IF (IRET.EQ.0) THEN IF (ILOCB.GE.0) THEN BASEN = VPARM(1+ILOCB) IA1 = BASEN / 256. + 0.1 IA2 = BASEN - IA1*256. + 0.1 ELSE IA1 = VPARM(1+ILOCA1) + 0.1 IA2 = VPARM(1+ILOCA2) + 0.1 END IF IF (.NOT.REQBAS (IA1, IA2, DESEL, IXANT, NXANT, IXBAS, NXBAS)) * GO TO 100 C Copy to output. XCOUNT = XCOUNT + 1.0D0 CALL RCOPY (NRPARM, VPARM, BUFF2(IPTRO)) CALL RCOPY (NCOPY, BUFF1, BUFF2(IPTRO+NRPARM)) IPTRO = IPTRO + LREC NIOUT = NIOUT + 1 C Write vis record. IF (NIOUT.GE.NIOLIM) THEN CALL UVDISK ('WRIT', LUN, FIND, BUFF2, NIOLIM, BIND, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'WRITING VIS DATA' GO TO 990 END IF IPTRO = BIND NIOUT = 0 END IF C Read next buffer. GO TO 100 END IF C Finish write NIOUT = - NIOUT CALL UVDISK ('FLSH', LUN, FIND, BUFF2, NIOUT, BIND, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'FINISHING WRITING VIS DATA' GO TO 990 END IF C tell NVIS = XCOUNT + 0.001D0 WRITE (MSGTXT,1100) NVIS IF (NVIS.LE.0) THEN IRET = 10 GO TO 990 ELSE CALL MSGWRT (5) IRET = 0 END IF C Compress output file. CALL UCMPRS (NVIS, WVOL, WCNO, LUN, CATBLK, IRET) C Close files CALL UVGET ('CLOS', VPARM, BUFF1, IRET) CALL ZCLOSE (LUN, FIND, IRET) C copy keywords CALL KEYCOP (DISKIN, OLDCNO, WVOL, WCNO, IERR) C Copy tables CALL COPTAB (DISKIN, OLDCNO, WVOL, WCNO, IERR) IF (IERR.GT.2) THEN MSGTXT = 'UVCOP: ERROR COPYING TABLES TO OUTPUT UV SCRATCH' CALL MSGWRT (6) END IF GO TO 999 C error 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('UVCOP: ERROR',I4,' ON ',A) 1100 FORMAT ('Copied',I12,' visibilities to the work file') END SUBROUTINE UVMDVS (APCORE, DISKI, CNOSCI, DISK1, CNOSC1, DISK2, * CNOSC2, MODEL, METHOD, DOMSG, CHANEL, NCHAN, JBUFSZ, FREQID, * BUFF1, BUFF2, BUFF3, BUFF4, BUFF5, IRET) C----------------------------------------------------------------------- C UVMDIV divides model visibilities derived from CLEAN components C or images into a uv data set. The weights of the data returned C will be the input values multiplied by the model amplitude. C A variety of model computation methods are available; if a single C pass thru UVMDVS, the DFT routine, is not sufficient then the data C is copied to a scratch file which has space for a second copy of the C data, the model values are computed and summed in these locations C and finally the model is divided into the data and written to the C output file. C Extensive use is made of commons to communicate with UVMDIV, in C particular /MAPDES/ (include DGDS.INC) contains most of the C critical information about the CLEAN components files or images to C be used. Common /UVHDR/ (filled in by UVPGET) is presumed to C describe the uv data files. C Also fills in frequency table (NCHANG, FREQG) in include C DGDS.INC C Inputs: C DISKI I Input disk number. if <= 0 then input is a C scratch file. C CNOSCI I Input file catalog slot number or /CFILES/ C scratch file number. C DISK1 I Output disk number. if <= 0 then output is a C scratch file. C CNOSC1 I Output file catalog slot number or /CFILES/ C scratch file number. If <= 0 then one of the C internal scratch files will be used. C MODEL I 1=> clean components, 2=>image. C METHOD I 1=>gridded, -1=>DFT, 0=>chose. C DOMSG L If true give percent done messages for DFT. C CHANEL I First uv data channel to subtract. C NCHAN I Number of frequency channels to subtract. C CATBLK I(256) UV data catalog header record. C JBUFSZ I Size of BUFF1,2,3 in bytes, must be at least C 4096 words. C FREQID I Freq ID number, if it exists. C Inputs from COMMON /MAPDES/: C MFIELD I Number of fields C NSUBG I(*) Number of components already sub. C NCLNG I(*) Number of components per field. C CCDISK I(*) Disk numbers for CC files C CCCNO I(*) Catalog slot numbers for CC files. C CCVER I(*) CC file version number for each field. C FACGRD R(2) Value to multiply clean component fluxes by C before subtraction (negative for sum). C FACGRD(2) is for data and 0 or 1 only values C used. Model added not subtracted when data are C ignored. C SCTYPE C*2 Scratch file type to create. (eg. 'SC') C NONEG L Stop reading comps. from a file past the first C negative component. C LIMFLX R Stop when flux < LIMFLX. C DOPTMD L Use the point model specified by PTFLX, PTRAOF, C PTDCOF (DFT modeling ONLY) C PTFLX R Point model flux density (Jy) (I pol. only) C PTRAOF R Point model RA offset from uv phase center C (asec) C PTDCOF R Point model Dec. offset from uv phase center C Input from COMMON /UVHDR/: C LREC I Length of visibility record. C NVIS I Number of visibility records. C NRPARM I "Random" parameters before data, can be used C to skip observed values when computing model. C Output: C CNOSCO I Output file catalog slot number or /CFILES/ C scratch file number. Value returned if not C specified in call. C BUFF1 R(*) Work buffers. C BUFF2 R(*) Work buffers. C BUFF3 R(*) Work buffers. C BUFF4 R(*) Work buffers. C BUFF5 R(*) Work buffers. C IRET I Return error code. 0=>OK, otherwise failed. C----------------------------------------------------------------------- DOUBLE PRECISION APCORE(*) INTEGER DISKI, CNOSCI, DISK1, CNOSC1, DISK2, CNOSC2, MODEL, * METHOD, CHANEL, NCHAN, JBUFSZ, FREQID, BUFF4(*), IRET LOGICAL DOMSG REAL BUFF1(*), BUFF2(*), BUFF3(*), BUFF5(*) C INTEGER I, LENMOD, ISCR2, DISKX, SAVNRP, SAVLRC, DISK, CNO, LUN, * LUN2, INMETH, OUMETH, CNOX, LBIF, LEIF, CATSCR(256), TEMP(256), * LENBU, MXCMP, XNCC, APSIZ, SCFRW, IV, OV, NEED, MSGSAV, KAP, * LCHANL, LCHAN LOGICAL DOSUM, DODFT, F, LTEMP REAL BUFSZ, TIMDFT, TIMFFT INCLUDE 'INCS:PUVD.INC' DOUBLE PRECISION SFOFF(MAXIF) INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DAPM.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DSCD.INC' INCLUDE 'INCS:DMOD.INC' INCLUDE 'INCS:DCAT.INC' SAVE LUN, LUN2 DATA LENMOD /7/, LUN/27/, LUN2/28/ DATA F /.FALSE./ DATA SFOFF /MAXIF*0.0D0/ C----------------------------------------------------------------------- IRET = 0 NGRDAT = F C Store CATBLK for later use CALL COPY (256, CATBLK, SCRCAT) SCLREC = LREC SCRPRM = NRPARM COMPDT = CATBLK(KINAX).EQ.1 IF (COMPDT) THEN CALL AXEFND (8, 'WEIGHT ', SCRCAT(KIPCN), SCRHOL(KHPTP), * WTLOC, IRET) IF ((IRET.NE.0) .OR. (WTLOC.LT.0)) THEN IRET = 5 MSGTXT = 'CANNOT FIND WEIGHT AND SCALE FOR COMPRESSED DATA' CALL MSGWRT (8) GO TO 999 END IF IRET = 0 END IF C Decide model computation method. INMETH = METHOD DODFT = (METHOD.LT.0) .AND. (MODEL.EQ.1) C DFT only for pt. model CALL UVMTYP (0, METHOD, NCHAN, LTEMP, TIMDFT, TIMFFT) IF (MODEL.EQ.1) DODFT = LTEMP DODFT = DODFT .OR. DOPTMD OUMETH = 1 IF (DODFT) OUMETH = -1 C If doing Direct Fourier Trans. IF (DODFT) THEN C Check point model IF (DOPTMD) THEN XNCC = 1 ELSE XNCC = 0 DO 20 I = 1,MFIELD XNCC = XNCC + NCLNG(I) - NSUBG(I) 20 CONTINUE END IF C Make AP memory for this BUFSZ = JBUFSZ LENBU = ((BUFSZ-2*NBPS) / 2) / (LREC*2) NEED = 12 + (LENBU*LREC) + CATBLK(KINAX+JLOCF) + XNCC*LENMOD * + 2*LENMOD NEED = NEED / 1024 MSGSAV = MSGSUP MSGSUP = 32000 CALL QINIT (APCORE, NEED, 0, KAP) MSGSUP = MSGSAV IF ((KAP.EQ.0) .OR. (PSAPNW.EQ.0)) THEN NEED = 12 + (LENBU*LREC) + CATBLK(KINAX+JLOCF) + * XNCC*LENMOD/10 + 2*LENMOD NEED = NEED / 1024 NEED = MIN (32+1024, NEED) + 2 CALL QINIT (APCORE, NEED, 0, KAP) IF ((KAP.EQ.0) .OR. (PSAPNW.EQ.0)) THEN IRET = 8 MSGTXT = 'UVMDIV CANNOT GET NEEDED MEMORY' CALL MSGWRT (8) END IF END IF APSIZ = PSAPNW * 1024 CALL QRLSE C Decide if will fit.in AP MXCMP = (APSIZ-10.-(LENBU*LREC)-2.-CATBLK(KINAX+JLOCF)) * / LENMOD MXCMP = MXCMP - 2 C If room for all components IF (XNCC.LE.MXCMP) THEN C IF using a scratch file IF (DISKI.LE.0) THEN DISK = SCRVOL(CNOSCI) CNO = SCRCNO(CNOSCI) ELSE DISK = DISKI CNO = CNOSCI END IF C Fill Frequency table. CALL FRQTAB (DISK, CNO, LUN, CATBLK, FREQID, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Divide model. CALL EVADFT (APCORE, CHANEL, NCHAN, DISKI, CNOSCI, DISK1, * CNOSC1, DISK2, CNOSC2, 0, DOSUM, DOMSG, CATR, JBUFSZ, * BUFF1,BUFF2, BUFF3, BUFF4, IRET) C Check for too many comps. IF (IRET.EQ.10) GO TO 100 GO TO 999 END IF END IF C Use UVMSUB for model. C Copy to padded scratch file. 100 ISCR2 = 0 C Message to about division MSGTXT = 'Divide data by model - first compute model by summing' IF (MSGSUP.NE.32000) CALL MSGWRT (3) CALL UVDPAD (DISKI, CNOSCI, ISCR2, JBUFSZ, BUFF1, BUFF2, IRET) IF (IRET.NE.0) GO TO 999 C Get scratch file header IF (DISKI.LE.0) THEN DISKX = SCRVOL(CNOSCI) CNOX = SCRCNO(CNOSCI) ELSE DISKX = DISKI CNOX = CNOSCI END IF CALL CATIO ('READ', DISKX, CNOX, CATSCR, 'REST', BUFF1, IRET) IF ((IRET.GE.1) .AND. (IRET.LE.4)) GO TO 999 IRET = 0 C Determine LBIF, LEIF LBIF = 1 LEIF = 1 C If more than 1 IF IF (JLOCIF.GT.0) THEN LEIF = CATBLK(KINAX+JLOCIF) C Copy part portion of IF table IV = 1 OV = 1 CALL CHNCOP (IV, OV, LUN, LUN2, DISKX, SCRVOL(ISCR2), CNOX, * SCRCNO(ISCR2), CATBLK, CATSCR, LBIF, LEIF, FREQID, SFOFF, * BUFF1, BUFF2, BUFF3, BUFF5, IRET) IF (IRET.NE.0) GO TO 999 C End if more than 1 IF END IF C Compute model. C Set factor for subtraction. FACGRD(1) = - FACGRD(1) C Redefine record size in C /UVHDR/. DISKX = 0 SAVNRP = NRPARM SAVLRC = LREC NRPARM = LREC LREC = 2 * SAVLRC - SAVNRP CALL UVMSUB (APCORE, DISKX, ISCR2, DISKX, ISCR2, 0, MODEL, METHOD, * CHANEL, NCHAN, F, DOMSG, CATBLK, JBUFSZ, FREQID, BUFF1, BUFF2, * BUFF3, IRET) IF (IRET.NE.0) GO TO 999 C Restore /UVHDR/ NRPARM = SAVNRP LREC = SAVLRC C Reset factor for subtraction. FACGRD(1) = - FACGRD(1) C Divide/compress record to output LCHANL = CHANEL LCHAN = NCHAN IF (MODMAX.GT.0) THEN LCHAN = 0 LCHANL = 1000000 DO 110 I = 1,MODMAX IV = MODCHN(I) + MODNCH(I) - 1 LCHAN = MAX (LCHAN, IV) LCHANL = MIN (LCHANL, MODCHN(I)) 110 CONTINUE LCHAN = LCHAN - LCHANL + 1 END IF CALL UVDSOU (ISCR2, DISK1, CNOSC1, DISK2, CNOSC2, LCHANL, * LCHAN, JBUFSZ, BUFF1, BUFF2, BUFF3, IRET) IF (IRET.NE.0) GO TO 999 C Delete scratch file IF (ISCR2.GT.0) THEN CALL COPY (256, CATBLK, TEMP) SCFRW = 2 IV = 1 CALL MAPCLR (IV, SCRVOL(ISCR2), SCRCNO(ISCR2), SCFRW, BUFF1) IF (ISCR2.EQ.NSCR) NSCR = NSCR - 1 CALL COPY (256, TEMP, CATBLK) END IF C 999 RETURN END SUBROUTINE UVDSOU (ISCR1, DISK1, CNOSC1, DISK2, CNOSC2, CHANEL, * NCHAN, JBUFSZ, BUFF1, BUFF2, BUFF3, IRET) C----------------------------------------------------------------------- C UVDSOU works on visibility data in which the lower half has C observed values and the upper half, in identical form, has model C values. The model values are divided into the observed values and C only the ratios are written to the output. C The weight of the output record is multiplied by the amplitude C of the model visibility. C The contents of /UVHDR/ are assumed valid for the output. C If disk .LE. 0 then CNOSCR is assumed to be the /CFILES/ scratch C file number. Uses LUNs 24 and 25. C Inputs: C ISCR1 I /CFILES/ number of scratch file for input C DISK1 I Disk number for catalogd output file. If .LE. 0 C then the output file is a /CFILES/ scratch file. C CNOSC1 I Catalog slot number if catalogd output file; C /CFILES/ scratch file number if a scratch file, C DISK2 I Disk number for catalogd output file. If .LE. 0 C then the output file is a /CFILES/ scratch file. C CNOSC2 I Catalog slot number if catalogd output file; C /CFILES/ scratch file number if a scratch file, C CHANEL I First uv data channel to divide. C NCHAN I Number of model channels to divide. C JBUFSZ I The size of BUFF1 and BUFF2 in bytes. C Inputs from common /UVHDR/: C LREC I Length of visibility record in words. C NRPARM I Number of random parameters C NVIS I Number of visibilities C Output: C BUFF1 R(*) Work buffers C BUFF2 R(*) Work buffers C ISCR2 I /CFILES/ number of output scratch file. C IRET I Return code, 0 => OK, otherwise failed. C----------------------------------------------------------------------- REAL BUFF1(*), BUFF2(*), BUFF3(*) INTEGER ISCR1, DISK1, CNOSC1, DISK2, CNOSC2, CHANEL, NCHAN, * JBUFSZ, IRET C INCLUDE 'INCS:ZPBUFSZ.INC' CHARACTER NAME*48 INTEGER VOL, FINDI, FIND1, FIND2, KLOOP, LUNI, LUN1, LUN2, * BINDI, BIND1, BIND2, LENBO, INIO, ONIO, LRECI, VVADD, OCNT, * ILOOP, INCVIS, JLOOP, LENBI, IPM, IPO, OP1, OP2, IIPM, IIPO, * OOP1, OOP2, BO, VO, LDATA, LRECUN, UNPARM LOGICAL T, F REAL AMP, AMPI, RATRE, RATIM, WT, TEMPBF(UVBFSS), WTSCL(2) INCLUDE 'INCS:PUVD.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DSCD.INC' DATA LUNI, LUN1, LUN2, BO, VO /24,25,26, 1,0/ DATA T, F /.TRUE.,.FALSE./ C----------------------------------------------------------------------- C Set output file name. IF (DISK1.LE.0) THEN VOL = SCRVOL(CNOSC1) CALL ZPHFIL ('SC', SCRVOL(CNOSC1), SCRCNO(CNOSC1), 1, NAME, * IRET) ELSE VOL = DISK1 CALL ZPHFIL ('UV', VOL, CNOSC1, 1, NAME, IRET) END IF C Open output file. CALL ZOPEN (LUN1, FIND1, VOL, NAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) 'WRIT', IRET GO TO 990 END IF C Set output file name. IF (DISK2.LE.0) THEN VOL = SCRVOL(CNOSC2) CALL ZPHFIL ('SC', SCRVOL(CNOSC2), SCRCNO(CNOSC2), 1, NAME, * IRET) ELSE VOL = DISK2 CALL ZPHFIL ('UV', VOL, CNOSC2, 1, NAME, IRET) END IF C Open output file. CALL ZOPEN (LUN2, FIND2, VOL, NAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) 'WRIT', IRET GO TO 990 END IF C Setup to init I/O, determine C size. LENBO = 0 CALL UVINIT ('WRIT', LUN1, FIND1, NVIS, VO, LREC, LENBO, JBUFSZ, * BUFF1, BO, BIND1, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1010) 'WRIT', IRET GO TO 990 END IF LENBO = 0 CALL UVINIT ('WRIT', LUN2, FIND2, NVIS, VO, LREC, LENBO, JBUFSZ, * BUFF2, BO, BIND2, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1010) 'WRIT', IRET GO TO 990 END IF OCNT = 0 ONIO = LENBO C Open input file. CALL ZPHFIL ('SC', SCRVOL(ISCR1), SCRCNO(ISCR1), 1, NAME, IRET) CALL ZOPEN (LUNI, FINDI, SCRVOL(ISCR1), NAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) 'READ', IRET GO TO 990 END IF C Setup for init. LRECI = 2 * SCLREC - SCRPRM IF (COMPDT) LRECI = SCLREC LENBI = 0 CALL UVINIT ('READ', LUNI, FINDI, NVIS, VO, LRECI, LENBI, JBUFSZ, * BUFF3, BO, BINDI, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1010) 'READ', IRET GO TO 990 END IF C Set up for reformatting. INCVIS = MIN (INCS, INCF) INCVIS = MAX (2, INCVIS) VVADD = (CHANEL - 1) * INCF + VOFF LRECUN = LREC UNPARM = NRPARM IF (COMPDT) THEN LRECUN = LREC - NRPARM LRECUN = LRECUN * 3 LRECUN = LRECUN + NRPARM - 2 UNPARM = NRPARM - 2 END IF LDATA = (LRECUN - UNPARM) / INCS IF (LDATA.GT.UVBFSS) THEN MSGTXT = 'UVDSOU: TEMPORARY BUFFER TOO SMALL - MODIFY' IRET = 1 GO TO 990 END IF C Begin loop 100 CALL UVDISK ('READ', LUNI, FINDI, BUFF3, INIO, BINDI, IRET) IF (INIO.LE.0) GO TO 250 IF (IRET.NE.0) THEN WRITE (MSGTXT,1100) 'READ', IRET GO TO 990 END IF C Loop copying to output, zero C padding. DO 200 ILOOP = 1,INIO C Copy random parameters and obs. CALL RCOPY (LRECUN, BUFF3(BINDI), BUFF1(BIND1)) CALL RCOPY (LRECUN, BUFF3(BINDI), BUFF2(BIND2)) IIPO = BINDI + UNPARM + VVADD IIPM = BINDI + LRECUN + VVADD OOP1 = BIND1 + UNPARM + VVADD OOP2 = BIND2 + UNPARM + VVADD C Divide C Loop over channel. DO 150 JLOOP = 1,NCHAN C Set indices. IPO = IIPO IPM = IIPM OP1 = OOP1 OP2 = OOP2 C Loop over Stokes DO 130 KLOOP = 1,NSTOK C subtraction BUFF1(OP1) = BUFF3(IPO) - BUFF3(IPM) BUFF1(OP1+1) = BUFF3(IPO+1) - BUFF3(IPM+1) BUFF1(OP1+2) = BUFF3(IPO+2) C division AMP = BUFF3(IPM)*BUFF3(IPM) + * BUFF3(IPM+1)*BUFF3(IPM+1) RATRE = 0.0 RATIM = 0.0 WT = 0.0 C Check for small values. IF (AMP.GE.1.0E-20) THEN AMPI = 1.0 / AMP RATRE = AMPI * ( BUFF3(IPO)*BUFF3(IPM) + * BUFF3(IPO+1)*BUFF3(IPM+1)) RATIM = AMPI * ( BUFF3(IPM)*BUFF3(IPO+1) - * BUFF3(IPO)*BUFF3(IPM+1)) WT = AMP * BUFF3(IPO+2) END IF C Fill in. BUFF2(OP2) = RATRE BUFF2(OP2+1) = RATIM BUFF2(OP2+2) = WT C Update pointers for Stokes IPO = IPO + INCS IPM = IPM + INCS OP1 = OP1 + INCS OP2 = OP2 + INCS 130 CONTINUE C Update pointers for freq. IIPO = IIPO + INCF IIPM = IIPM + INCF OOP1 = OOP1 + INCF OOP2 = OOP2 + INCF 150 CONTINUE C Compress data if needed IF (COMPDT) THEN CALL ZUVPAK (LDATA, BUFF1(BIND1+UNPARM), * WTSCL, TEMPBF) CALL RCOPY (2, WTSCL, BUFF1(BIND1+UNPARM)) CALL RCOPY (LDATA, TEMPBF, BUFF1(BIND1+UNPARM+2)) CALL ZUVPAK (LDATA, BUFF2(BIND2+UNPARM), * WTSCL, TEMPBF) CALL RCOPY (2, WTSCL, BUFF2(BIND2+UNPARM)) CALL RCOPY (LDATA, TEMPBF, BUFF2(BIND2+UNPARM+2)) END IF C Update pointers. BINDI = BINDI + LRECI BIND1 = BIND1 + LREC BIND2 = BIND2 + LREC C Check if time for output. OCNT = OCNT + 1 C Write IF (OCNT.GE.ONIO) THEN ONIO = OCNT CALL UVDISK ('WRIT', LUN1, FIND1, BUFF1, ONIO, BIND1, * IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1100) 'WRIT', IRET GO TO 990 END IF ONIO = OCNT CALL UVDISK ('WRIT', LUN2, FIND2, BUFF2, ONIO, BIND2, * IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1100) 'WRIT', IRET GO TO 990 END IF OCNT = 0 END IF 200 CONTINUE C Loop back for more. GO TO 100 C Done - flush output buffer. 250 ONIO = -OCNT CALL UVDISK ('FLSH', LUN1, FIND1, BUFF1, ONIO, BIND1, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1100) 'FLSH', IRET GO TO 990 END IF ONIO = -OCNT CALL UVDISK ('FLSH', LUN2, FIND2, BUFF2, ONIO, BIND2, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1100) 'FLSH', IRET GO TO 990 END IF C Close files. CALL ZCLOSE (LUNI, FINDI, IRET) CALL ZCLOSE (LUN1, FIND1, IRET) CALL ZCLOSE (LUN2, FIND2, IRET) IRET = 0 GO TO 999 C Error 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('UVDSOU: OPEN FOR ',A4,' ERROR ',I5) 1010 FORMAT ('UVDSOU: INIT FOR ',A4,' ERROR ',I5) 1100 FORMAT ('UVDSOU: ',A4,' ERROR ',I5) END SUBROUTINE EVADFT (APCORE, CHANEL, NCHAN, DISKI, CNOSCI, DISK1, * CNOSC1,DISK2, CNOSC2, IFIELD, DOSUM, DOMSG, CATR, JBUFSZ, * BUFF1, BUFF2, BUFF3, IBUFF, IRET) C----------------------------------------------------------------------- C EVADFT subtracts AND divides CLEAN components from/into ungridded C visibility data by a DFT model computation. Only model components C of a single type are processed. Point components will be taken as C Gaussians or Spheres as needed if some of the fields are extended C and some not. C All un subtracted data processed in one call. C Inputs: C CHANEL I Frequency channel: used if MODMAX=0 C NCHAN I Number of frequency channels.: ditto C DISKI I Input file disk number for cataloged files, C .LE. 0 => /CFILES/ scratch file. C CNOSCI I Input file catalog slot number or /CFILES/ C scratch file number. C DISK1 I Output file disk number for cataloged files, C .LE. 0 => /CFILES/ scratch file. subtraction C DISK2 I Output file disk number for cataloged files, C .LE. 0 => /CFILES/ scratch file. division C IFIELD I Field to do (0 -> all): used if MODMAX=0 C DOSUM L If true sum the flux in each field C DOMSG L If true give percent done messages. C CATR R(256) UV data catalog header record. C JBUFSZ I Size of BUFF1,2, IBUFF in AIPS bytes, each C must be at least 4096 words. C Inputs: from commons C MODMAX I DMOD.INC - if set, this controls facets/chans C MFIELD I Number of fields C NCLNG I(16) Number of components per field. - C changed if flux limit hit C NSUBG I(16) The next component to subtract. C CCDISK I(16) Disk numbers of the clean images. C CCCNO I(16) Catalog slot numbers of clean images. C CCVER I(*) CC file version number for each field. C NGRDAT L If FALSE get map size, scaling etc. parms C from the model map cat. header. If TRUE C then the values filled in by GRDAT must C already be filled into the common. C LREC I Length in words of vis record. C NVIS I Number of vis. records C NONEG L Stop reading comps. from a file past the first C negative component. C LIMFLX R Stop if abs(flux) < LIMFLX C DOPTMD L Use the point model specified by PTFLX, PTRAOF, C PTDCOF C PTFLX R Point model flux density (Jy) (I pol. only) C PTRAOF R Point model RA offset from uv phase center C (asec) C PTDCOF R Point model Dec. offset from uv phase center C PARMOD R(6) Model parameters for non point models; used C only if DOPTMOD is true. C 1=> model type, 0=point, 1=gaussian, 3=sphere C Gaussian: (2)=major axis(asec), (3)=minor axis C (4)=PA (degrees) C Sphere: (2)=radius (asec). C KSTOK I (DGDS.INC) If a point model is specified a C value of 2 indicates a Q pol model, 3 U and C 4 V pol.AC C In/out: C CNOSC1 I IN: output file catalog slot number or /CFILES/ C scratch file number. Will create a scratch file C if CNOSC1 and DISK1 .le. 0. C Out: file /CFILES/ number if created. C CNOSC2 I Output file 2 catalog slot number or /CFILES/ C scratch file number. C Output: C BUFF1 R(*) I/O buffers. C BUFF2 R(*) I/O buffers. C IBUFF I(*) I/O Buffer. C IRET I Return code, 0 => ok, otherwise not. C 8 => model types not compatible/illegal C 9 => Buffers too small to load AP. C 10 => Too many components for division. C----------------------------------------------------------------------- DOUBLE PRECISION APCORE(*) INTEGER CHANEL, NCHAN, DISKI, CNOSCI, DISK1, CNOSC1, DISK2, * CNOSC2, IFIELD, JBUFSZ, IBUFF(*), IRET LOGICAL DOSUM, DOMSG REAL BUFF1(*), BUFF2(*), BUFF3(*), CATR(256) C CHARACTER KEYS(7)*8, NAME*48, MDTYP(4)*8, ERRTXT*40, UMET*4 INTEGER JCOMP, JNCOMP, CCOUNT, XNCOMP, MXCMP, CURCMP, JT, * NCOMP, J, MCOMP, VO, BO, ISIZE, INIO, MMCOMP, NNCOR, IDATA, UV, * LLREC, IAPBUF, IAPCC0, IAPCT, LMCOMP, IAPTMP, VIS, WRK, LLNMOD, * MCHAN, JNCS, JNCF, KAP, SFLAG, APSIZ, MXCC, INIO2, LUNC, VOL, * INDEX, ITYPE, NIOUT, KBIND1, KBIND2, LENBU, LENMOD, JLREC, * JNREC, FINDI, FIND1, FIND2, I, LUNI, LUN1, LUN2, ITIME(3), * IBIND, LFIELD, LMOD(4), NKEY, KOLS(7), RAKOL, DECKOL, FLXKOL, * TYPKOL, MAJKOL, MINKOL, PAKOL, IPCLST, IPCDNE, NTIMES, MODTYP, * LRPARM, LF1, LF2, NEED, MSGSAV, LCHAN, I1, I2, FLD1, FLD2, * CHN1, CHN2 REAL XXOFF, YYOFF, ZZOFF, RECORD(20), FACT2(4), CPA, SPA, * XMAJ, XMIN, ABFACG, XYZ(3), XP(3), UMAT(3,3), PMAT(3,3), UUU, * BTEMP DOUBLE PRECISION XTLST, PCTOT, PCLST, XRA, XDEC, XPR, YPR, CONST, * CONST2 LOGICAL T, F, WESET INCLUDE 'INCS:PSTD.INC' INCLUDE 'INCS:PUVD.INC' INCLUDE 'INCS:DMOD.INC' LOGICAL DONE(MAXFLD), ONZE INTEGER IAPCC, MFR, NFR, LFR, MMCMP(MAXMOD) DOUBLE PRECISION FFRAC, X, Y, PTFL0 REAL FBUFF(MAXCIF) INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DMPR.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DAPM.INC' INCLUDE 'DSPX.INC' PARAMETER (CONST = DG2RAD * TWOPI) SAVE ONZE EQUIVALENCE (KOLS(1), RAKOL), (KOLS(2), DECKOL), * (KOLS(3), FLXKOL), (KOLS(4), TYPKOL), * (KOLS(5) ,MAJKOL), (KOLS(6), MINKOL), * (KOLS(7), PAKOL) DATA KEYS /'DELTAX ', 'DELTAY ', 'FLUX ', * 'TYPE OBJ', 'MAJOR AX', 'MINOR AX', 'POSANGLE'/ DATA LMOD /4, 7, 7, 6/ DATA MDTYP /'Point ', 'Gaussian', 'Unknown ', 'Sphere '/ DATA VO, BO, MXCC /0, 1, 1024/ DATA LUNI, LUN1, LUN2, LUNC /22,23,24,29/ DATA T, F /.TRUE.,.FALSE./, ONZE/.FALSE./ C----------------------------------------------------------------------- UMET = 'DFT' CALL FILL (MAXMOD, 0, MMCMP) C CONST2 converts FWHM(deg) to C coefficients for u*u, v*v, u*v CONST2 = DG2RAD * (PI / 1.17741022D0) * SQRT (0.5D0) C Tell kind of operation. MSGTXT = 'EVADFT: Begin DFT component subtraction & division' CALL MSGWRT (2) NNCOR = 1 C Get un-compressed UV increments CALL UVINCS (INCS, INCF, INCIF, NRPARM, LREC, JNCS, JNCF, LRPARM, * LLREC) LFIELD = 0 CCOUNT = 0 IF (IFIELD.LE.0) THEN LF1 = 1 LF2 = MFIELD ELSE LF1 = IFIELD LF2 = IFIELD END IF WESET = (MODMAX.LE.0) .AND. (.NOT.DOPTMD) MCHAN = NCHAN IF (WESET) THEN MODMAX = LF2 - LF1 + 1 DO 5 I = 1,MODMAX MODFLD(I) = I + LF1 - 1 MODCHN(I) = CHANEL MODNCH(I) = NCHAN MODCCV(I) = ABS(CCVER(MODFLD(I))) MODCCB(I) = NSUBG(MODFLD(I)) 5 CONTINUE FLD1 = LF1 FLD2 = LF2 CHN1 = CHANEL CHN2 = CHANEL + NCHAN - 1 ELSE IF (MODMAX.GT.0) THEN MCHAN = 0 FLD1 = 100000 FLD2 = 0 CHN1 = FLD1 CHN2 = 0 DO 6 I = 1,MODMAX MCHAN = MAX (MCHAN, MODNCH(I)) FLD1 = MIN (FLD1, MODFLD(I)) FLD2 = MAX (FLD2, MODFLD(I)) CHN1 = MIN (CHN1, MODCHN(I)) CHN2 = MAX (CHN2, MODCHN(I)+MODNCH(1)-1) 6 CONTINUE END IF IF (MODMAX.GT.0) THEN WRITE (MSGTXT,1006) FLD1, FLD2, CHN1, CHN2, MODMAX CALL REFRMT (MSGTXT, ' ', I) CALL MSGWRT (2) END IF C Decide component type. C From model passed IF (DOPTMD) THEN MODTYP = PARMOD(1) + 0.5 IF (MODTYP.EQ.2) MODTYP = 3 XNCOMP = 1 MODMAX = 1 MODFLD(1) = LF1 MODCHN(1) = CHANEL MODNCH(1) = NCHAN MODCCV(1) = 0 MODCCB(1) = 1 MMCMP(1) = 1 WESET = .TRUE. IF (DOPTSP) MCHAN = 1 C From CC table, field 1 ELSE LFIELD = LF1 C Get field info. if nec. IF (.NOT.NGRDAT) THEN CALL GRDAT (F, LFIELD, CATR, IBUFF(2049), IRET) IF (IRET.NE.0) GO TO 999 END IF C If NGRDAT read CLEAN CATBLK. IF (NGRDAT) THEN ERRTXT = 'READING CLEAN CATBLK' CALL CATIO ('READ', CCDISK(LFIELD), CCCNO(LFIELD), KLNBLK, * 'REST', IBUFF(2049), IRET) IF ((IRET.GT.0) .AND. (IRET.LT.5)) GO TO 990 END IF C For point model MODTYP = 0 XNCOMP = 0 C check all CC files MFR = 0 LFR = 0 DO 10 I = 1,MODMAX IF (MODCHN(I).NE.LFR) THEN MFR = MFR + 1 LFR = MODCHN(I) END IF LFIELD = MODFLD(I) JNREC = 1 JLREC = 0 NKEY = 0 ERRTXT = 'OPENING CLEAN COMPS FILE' CALL TABINI ('READ', 'CC', CCDISK(LFIELD), CCCNO(LFIELD), * MODCCV(I), KLNBLK, LUNC, NKEY, JNREC, JLREC, BUFF1, * IBUFF, IRET) IF (IRET.GT.1) GO TO 990 IF (NCLNG(LFIELD).LE.0) NCLNG(LFIELD) = IBUFF(5) DONE(I) = F XNCOMP = XNCOMP + NCLNG(LFIELD) - MODCCB(I) + 1 MMCMP(MFR) = MMCMP(MFR) + NCLNG(LFIELD) - MODCCB(I) + 1 C More complex models C Find columns (physical) IF ((JLREC.GT.3) .AND. (NCLNG(LFIELD).GE.MODCCB(I))) * THEN NKEY = JLREC IF (JLREC.GT.7) NKEY = 7 ERRTXT = 'FINDING REQUIRED CC COLUMNS' CALL FNDCOL (NKEY, KEYS, 8, F, IBUFF, KOLS, IRET) IF (IRET.NE.0) GO TO 990 C Read 1st record JCOMP = MODCCB(I) CALL TABIO ('READ', 0, JCOMP, RECORD, IBUFF, IRET) IF (IRET.GT.0) THEN WRITE (MSGTXT,1090) IRET, JCOMP GO TO 995 END IF C Get model type. ITYPE = RECORD(TYPKOL) + 0.5 IF (LFIELD.EQ.LF1) MODTYP = ITYPE IF (MODTYP.NE.ITYPE) THEN IF ((MODTYP.EQ.0) .OR. (ITYPE.EQ.0)) THEN MSGTXT = 'EVADFT: will treat points as extended' CALL MSGWRT(6) MODTYP = MAX (MODTYP, ITYPE) ELSE WRITE (MSGTXT,1000) MODTYP, ITYPE IRET = 8 GO TO 995 END IF END IF END IF C Close CLEAN components file. CALL TABIO ('CLOS', 0, JCOMP, RECORD, IBUFF, IRET) 10 CONTINUE END IF C Bad model type. IF ((MODTYP.NE.0) .AND. (MODTYP.NE.1) .AND. (MODTYP.NE.3)) THEN IRET = 8 WRITE (MSGTXT,1002) MODTYP GO TO 995 END IF C Check for point model. LFIELD = LF1 - 1 C Tell model type once IF (.NOT.ONZE) THEN C Tell model type MSGTXT = 'EVADFT: Model components of type '//MDTYP(MODTYP+1) CALL MSGWRT (2) ONZE = .NOT. ONZE C Check Buffer size IF ((XNCOMP.GT.10) .AND. (JBUFSZ/2.LT.4096)) THEN MSGTXT = 'EVADFT: SCRATCH BUFFER TOO SMALL FOR CCs!' CALL MSGWRT (8) IRET = 9 GO TO 999 END IF END IF C Set model length LENMOD = LMOD(MODTYP+1) C Determine size of uv I/O and C the number of CC that will fit. LENBU = ((JBUFSZ-2*NBPS) / 2) / (LLREC*2) C How much data fits in AP? JT = 15 + LENBU * LLREC IF (MXCC.GT.JT) JT = MXCC NEED = JT + MCHAN + (XNCOMP+2)*LENMOD NEED = NEED / 1024 MSGSAV = MSGSUP MSGSUP = 32000 CALL QINIT (APCORE, NEED, 0, KAP) MSGSUP = MSGSAV IF ((KAP.EQ.0) .OR. (PSAPNW.EQ.0)) THEN NEED = JT + MCHAN + (XNCOMP/10+2)*LENMOD NEED = NEED / 1024 NEED = MIN (32*1024, NEED) + 2 CALL QINIT (APCORE, NEED, 0, KAP) IF ((KAP.EQ.0) .OR. (PSAPNW.EQ.0)) THEN IRET = 8 MSGTXT = 'EVADFT CANNOT GET NEEDED MEMORY' GO TO 995 END IF END IF APSIZ = PSAPNW * 1024 CALL QRLSE MXCMP = (APSIZ - JT - NCHAN) / LENMOD MXCMP = MXCMP - 2 CALL RFILL (20, 0.0, RECORD) C Set AP pointers. C UV=UV pointer, VIS=vis pointer C IAPCC0=CLEAN components pointer. IDATA = 12 + MCHAN UV = IDATA + ILOCU IAPCC0 = 15 + LENBU * LLREC + NCHAN IAPCC0 = MAX (IAPCC0, MXCC) LLNMOD = LENMOD C Compute number of passes. NTIMES = (1.0 * XNCOMP) / MXCMP + 0.99999 NTIMES = MAX (NTIMES, 1) C Only one pass allowed for C division. No can do. IF (NTIMES.GT.1) THEN ERRTXT = 'TOO MANY COMPONENTS FOR DIVISION' IRET = 10 GO TO 990 END IF C Fix for Division scaling ABFACG = ABS(FACGRD(1)) C Open uv files. C Set input file name. IF (DISKI.LE.0) THEN VOL = SCRVOL(CNOSCI) CALL ZPHFIL ('SC', VOL, SCRCNO(CNOSCI), 1, NAME, IRET) ELSE VOL = DISKI CALL ZPHFIL ('UV', VOL, CNOSCI, 1, NAME, IRET) END IF C Open input file. ERRTXT = 'OPEN-FOR-READ VIS FILE' CALL ZOPEN (LUNI, FINDI, VOL, NAME, T, F, T, IRET) IF (IRET.NE.0) GO TO 990 C Create scratch file if necessary IF ((DISK1.LE.0) .AND. (CNOSC1.EQ.0)) THEN CALL UVSIZE (LREC, NVIS, ISIZE) ERRTXT = 'CREATING SCRATCH 1 FILE' CALL SCREAT (ISIZE, BUFF2, IRET) CNOSC1 = NSCR IF (IRET.GT.0) THEN IF (IRET.EQ.1) ERRTXT = 'NO SPACE FOR SCRATCH 1 FILE' GO TO 990 END IF END IF IF ((DISK2.LE.0) .AND. (CNOSC2.EQ.0)) THEN CALL UVSIZE (LREC, NVIS, ISIZE) ERRTXT = 'CREATING SCRATCH 2 FILE' CALL SCREAT (ISIZE, BUFF2, IRET) CNOSC2 = NSCR IF (IRET.GT.0) THEN IF (IRET.EQ.1) ERRTXT = 'NO SPACE FOR SCRATCH 2 FILE' GO TO 990 END IF END IF C Open vis file for write. IF (DISK1.LE.0) THEN VOL = SCRVOL(CNOSC1) CALL ZPHFIL ('SC', VOL, SCRCNO(CNOSC1), 1, NAME, IRET) ELSE VOL = DISK1 CALL ZPHFIL ('UV', VOL, CNOSC1, 1, NAME, IRET) END IF ERRTXT = 'OPEN-FOR-WRITE VIS 1 FILE' CALL ZOPEN (LUN1, FIND1, VOL, NAME, T, F, T, IRET) IF (IRET.NE.0) GO TO 990 IF (DISK2.LE.0) THEN VOL = SCRVOL(CNOSC2) CALL ZPHFIL ('SC', VOL, SCRCNO(CNOSC2), 1, NAME, IRET) ELSE VOL = DISK2 CALL ZPHFIL ('UV', VOL, CNOSC2, 1, NAME, IRET) END IF ERRTXT = 'OPEN-FOR-WRITE VIS 2 FILE' CALL ZOPEN (LUN2, FIND2, VOL, NAME, T, F, T, IRET) IF (IRET.NE.0) GO TO 990 C Loop, subtracting Max Component C In AP each pass I1 = 1 I2 = I1 - 1 C Setup for % done messages. PCTOT = NVIS PCLST = 0 IPCLST = 0 C Set AP loc for next CC load IAPCT = IAPCC0 C Determine. no. this pass. MMCOMP = 0 NCOMP = MIN (MXCMP, XNCOMP) C Grab AP. CALL QINIT (APCORE, NEED, 0, KAP) IF ((KAP.EQ.0) .OR. (PSAPNW.LE.0)) THEN IRET = 10 MSGTXT = 'EVADFT: BIZARRE FAILURE TO GET AP MEMORY' GO TO 995 END IF C Initialize REAL time clock for C AP roller. CALL ZTIME (ITIME) XTLST = (ITIME(1) * 60.00) + ITIME(2) + (ITIME(3) / 60.0) C If Not single component model. IF (.NOT.DOPTMD) THEN C Find next FIELD. C Loop back here for next field. 70 I2 = I2 + 1 C See if done. IF (I2.GT.MODMAX) GO TO 150 IF (DONE(I2)) GO TO 70 LFIELD = MODFLD(I2) C See if there are CCs. IF ((MODCCB(I2).GT.NCLNG(LFIELD))) GO TO 70 C Get field info. if nec. IF (.NOT.NGRDAT) THEN CALL GRDAT (F, LFIELD, CATR, IBUFF(2049), IRET) IF (IRET.NE.0) GO TO 999 END IF C If NGRDAT read CLEAN CATBLK. IF (NGRDAT) THEN ERRTXT = 'READING CLEAN CATBLK' CALL CATIO ('READ', CCDISK(LFIELD), CCCNO(LFIELD), * KLNBLK, 'REST', IBUFF(2049), IRET) IF ((IRET.GT.0) .AND. (IRET.LT.5)) GO TO 990 END IF C Set field center offsets. XXOFF = DXCG(LFIELD) * CCROT + DYCG(LFIELD) * SSROT YYOFF = DYCG(LFIELD) * CCROT - DXCG(LFIELD) * SSROT ZZOFF = DZCG(LFIELD) CALL XYSHFT (RA, DEC, XSHIFT(LFIELD), YSHIFT(LFIELD), * MAPROT, XRA, XDEC) IF (DO3DIM) THEN CALL PRJMAT (RA, DEC, UVROT, XRA, XDEC, MAPROT, UMAT, * PMAT) ELSE CALL P2DMAT (RA, DEC, UVROT, XRA, XDEC, MAPROT, UMAT, * PMAT) END IF C Load CLEAN components into AP. C Open components file. JNREC = 1 JLREC = 0 NKEY = 0 ERRTXT = 'OPENING CLEAN COMPS FILE' CALL TABINI ('READ', 'CC', CCDISK(LFIELD), CCCNO(LFIELD), * MODCCV(I2), KLNBLK, LUNC, NKEY, JNREC, JLREC, BUFF1, IBUFF, * IRET) IF (IRET.GT.1) GO TO 990 JCOMP = MODCCB(I2) C Find columns (physical) NKEY = JLREC TYPKOL = 7 IF (JLREC.GT.7) NKEY = 7 ERRTXT = 'FINDING REQUIRED CC COLUMNS' CALL FNDCOL (NKEY, KEYS, 8, F, IBUFF, KOLS, IRET) IF (IRET.GT.0) GO TO 990 C Make sure that there are some IF (IBUFF(5).LE.0) GO TO 140 IF (NCLNG(LFIELD).LE.0) NCLNG(LFIELD) = IBUFF(5) C Loop loading components. IAPBUF = 10 CURCMP = MMCOMP + 1 C Check next component IF (JCOMP.GT.NCLNG(LFIELD)) GO TO 140 DO 130 J = CURCMP,NCOMP,MXCC JT = J - 1 JNCOMP = 0 MCOMP = NCOMP - J + 1 IF (MCOMP.GT.MXCC) MCOMP = MXCC IF (MCOMP.GT.(NCLNG(LFIELD)-MODCCB(I2)+1)) * MCOMP = NCLNG(LFIELD) - MODCCB(I2) + 1 DO 110 I = 1,MCOMP C Check if finished field IF (JCOMP.GT.NCLNG(LFIELD)) GO TO 120 CALL TABIO ('READ', 0, JCOMP, RECORD, IBUFF, IRET) IF (IRET.GT.0) THEN WRITE (MSGTXT,1090) IRET, JCOMP GO TO 995 END IF C Check that point comp. JCOMP = JCOMP + 1 JT = JT + 1 IF (JLREC.LE.3) THEN ITYPE = 0 ELSE ITYPE = RECORD(TYPKOL) + 0.5 END IF IF ((ITYPE.NE.MODTYP) .AND. (ITYPE.NE.0)) THEN WRITE (MSGTXT,1070) LFIELD, JCOMP-1 CALL MSGWRT (6) END IF IF (((ITYPE.EQ.MODTYP) .OR. (ITYPE.EQ.0)) .AND. * (IRET.EQ.0)) THEN C Check negative component limit DONE(I2) = (NONEG.AND.(RECORD(FLXKOL).LE.0.0)) * .OR. (ABS(RECORD(FLXKOL)).LT.LIMFLX) IF (DONE(I2)) THEN NCLNG(LFIELD) = JCOMP - 1 GO TO 120 END IF C If req. sum flux IF (DOSUM) THEN FLUXG(LFIELD) = FLUXG(LFIELD) + RECORD(FLXKOL) TFLUXG = TFLUXG + RECORD(FLXKOL) END IF JNCOMP = JNCOMP + 1 CCOUNT = CCOUNT + 1 XP(1) = (RECORD(RAKOL) + XPOFF(LFIELD)) * CONST XP(2) = (RECORD(DECKOL) + YPOFF(LFIELD)) * CONST XP(3) = 0.0 CALL PRJMUL (2, XP, UMAT, XYZ) BUFF1(JNCOMP) = XYZ(1) + XXOFF BUFF1(1024+JNCOMP) = XYZ(2) + YYOFF BUFF1(2048+JNCOMP) = XYZ(3) + ZZOFF C Handle scaling for division BUFF1(3072+JNCOMP) = ABFACG * RECORD(FLXKOL) C Gaussian IF (MODTYP.EQ.1) THEN C Convert to convenient C coefficients. IF (ITYPE.EQ.0) THEN BUFF2(JNCOMP) = 0.0 BUFF2(1024+JNCOMP) = 0.0 BUFF2(2048+JNCOMP) = 0.0 ELSE CPA = COS (DG2RAD*RECORD(PAKOL)) SPA = SIN (DG2RAD*RECORD(PAKOL)) XMAJ = RECORD(MAJKOL) * CONST2 XMIN = RECORD(MINKOL) * CONST2 BUFF2(JNCOMP) = - (((CPA * XMAJ)**2) + * (SPA * XMIN)**2) BUFF2(1024+JNCOMP) = - (((SPA * XMAJ)**2) * + (CPA * XMIN)**2) BUFF2(2048+JNCOMP) = - 2.0 * CPA * SPA * * (XMAJ*XMAJ - XMIN*XMIN) END IF C Sphere ELSE IF (MODTYP.EQ.3) THEN BUFF1(3072+JNCOMP) = 3.0 * BUFF1(3072+JNCOMP) IF (ITYPE.EQ.0) THEN BUFF2(JNCOMP) = 0.0 ELSE BUFF2(JNCOMP) = RECORD(MAJKOL) * 0.109662271 END IF BUFF2(1024+JNCOMP) = 0.1 END IF END IF 110 CONTINUE C Load components 120 IF (JNCOMP.GT.0) THEN LMCOMP = JNCOMP MMCOMP = MMCOMP + LMCOMP C Load into AP IAPBUF = 10 C x component CALL QPUT (APCORE, BUFF1, IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 1 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C y component CALL QPUT (APCORE, BUFF1(1025), IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 2 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C z component CALL QPUT (APCORE, BUFF1(2049), IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 3 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C Flux density CALL QPUT (APCORE, BUFF1(3073), IAPBUF, LMCOMP, 2) CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPCT, LLNMOD, LMCOMP) CALL QWR C Gaussian IF (MODTYP.EQ.1) THEN C Coef 1. CALL QPUT (APCORE, BUFF2, IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 4 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C Coef 2. CALL QPUT (APCORE, BUFF2(1025), IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 5 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C Coef 3. CALL QPUT (APCORE, BUFF2(2049), IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 6 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR END IF C Sphere IF (MODTYP.EQ.3) THEN C Radius CALL QPUT (APCORE, BUFF2, IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 4 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR C Minimum argument CALL QPUT (APCORE, BUFF2(1025), IAPBUF, LMCOMP, 2) IAPTMP = IAPCT + 5 CALL QWD CALL QVMOV (APCORE, IAPBUF, 1, IAPTMP, LLNMOD, LMCOMP) CALL QWR END IF IAPCT = IAPCT + (LLNMOD * LMCOMP) END IF C Check if finished field. IF ((JCOMP.GT.NCLNG(LFIELD)) .OR. DONE(I2)) GO TO 140 130 CONTINUE C Close CLEAN components file. 140 CALL TABIO ('CLOS', 0, JCOMP, RECORD, IBUFF, IRET) C Update field sub. count. NSUBG(LFIELD) = JCOMP C Check if need another field. IF (JT.LT.NCOMP) GO TO 70 C Check no. comps. 150 IF (MMCOMP.EQ.0) THEN ERRTXT = 'NO POINT COMPONENTS FOUND' IRET = 1 GO TO 990 END IF C Load correlator factors CALL GETCTL (CATR, FACT2, IRET) IF (IRET.NE.0) GO TO 999 C Else, single component model. ELSE I2 = 1 CALL RFILL (LENMOD, 0.0, BUFF1) BUFF1(2) = PTRAOF * CONST / 3600.D0 BUFF1(3) = PTDCOF * CONST / 3600.D0 C single comp uses W term XPR = PTRAOF / 206264.81D0 YPR = PTDCOF / 206264.81D0 BUFF1(4) = -(SQRT (1.0D0 - XPR*XPR - YPR*YPR) - 1.0D0) * * 206264.81D0 * CONST / 3600.D0 C Point IF (MODTYP.EQ.0) BUFF1(1) = PTFLX C Gaussian IF (MODTYP.EQ.1) THEN C Convert to convenient C coefficients. CPA = COS (DG2RAD * PARMOD(4)) SPA = SIN (DG2RAD * PARMOD(4)) XMAJ = PARMOD(2) * CONST2 * 2.77777778E-4 XMIN = PARMOD(3) * CONST2 * 2.77777778E-4 BUFF1(5) = -(((CPA * XMAJ)**2) + (SPA * XMIN)**2) BUFF1(6) = -(((SPA * XMAJ)**2) + (CPA * XMIN)**2) BUFF1(7) = -2.0 * CPA * SPA * (XMAJ*XMAJ - XMIN*XMIN) BUFF1(1) = PTFLX END IF C Uniform sphere IF (MODTYP.EQ.3) THEN BUFF1(5) = PARMOD(2) * 0.109662271 * 2.7777778E-4 BUFF1(1) = PTFLX * 3.0 BUFF1(6) = 0.1 END IF PTFL0 = BUFF1(1) MMCOMP = 1 IAPCC = IAPCC0 CALL QPUT (APCORE, BUFF1, IAPCC, LLNMOD, 2) IAPCT = IAPCC + LLNMOD C Set Stokes for point model FACT2(1) = 1.0 FACT2(2) = 1.0 IF (ICOR0.LT.0) THEN C RR,LL etc. C Q? IF (KSTOK.EQ.2) THEN FACT2(1) = 1.0 FACT2(2) = 1.0 VOFF = (3 - ABS (ICOR0)) * INCS ELSE IF (KSTOK.EQ.3) THEN FACT2(1) = 1.0 FACT2(2) = -1.0 VOFF = (3 - ABS (ICOR0)) * INCS ELSE IF (KSTOK.EQ.4) THEN FACT2(1) = 1.0 FACT2(2) = -1.0 VOFF = 0 END IF ELSE C True Stokes IF ((KSTOK.GE.2) .AND. (KSTOK.LE.4)) THEN FACT2(1) = 1.0 FACT2(2) = 0.0 VOFF = (KSTOK - ICOR0) * INCS END IF END IF NSTOK = 1 IF ((CATR(KRCIC+JLOCS).LT.0.0) .AND. (NCOR.GE.2)) NSTOK = 2 IF ((NSTOK.EQ.2) .AND. (ABS (ICOR0).EQ.2)) NSTOK = 1 C End if not single comp. model END IF C Correct for FACGRD NNCOR = NSTOK FACT2(1) = FACT2(1) * FACGRD(1) FACT2(2) = FACT2(2) * FACGRD(1) WRK = MCHAN + 2 UUU = -1.0 IF (FACGRD(2).EQ.0.0) UUU = 0.0 CALL QPUT (APCORE, UUU, WRK+NNCOR, 1, 2) CALL QPUT (APCORE, FACT2, WRK, NNCOR, 2) SFLAG = 1 C Set flag for UPOL and RL,LR data IF ((KSTOK.EQ.3) .AND. (ICOR0.LT.0)) SFLAG = -1 C Set vis pointer VIS = IDATA + LRPARM + (CHANEL-1) * JNCF + (VOFF/INCS)*JNCS CALL QWAIT C Init for read & write C visibility file C Init vis file for write ERRTXT = 'INIT-FOR-READ VIS FILE' CALL UVINIT ('READ', LUNI, FINDI, NVIS, VO, LREC, LENBU, JBUFSZ, * BUFF3, BO, IBIND, IRET) IF (IRET.NE.0) GO TO 990 C Init vis file for read. ERRTXT = 'INIT-FOR-WRITE VIS 1 FILE' CALL UVINIT ('WRIT', LUN1, FIND1, NVIS, VO, LREC, LENBU, JBUFSZ, * BUFF1, BO, KBIND1, IRET) IF (IRET.NE.0) GO TO 990 ERRTXT = 'INIT-FOR-WRITE VIS 2 FILE' CALL UVINIT ('WRIT', LUN2, FIND2, NVIS, VO, LREC, LENBU, JBUFSZ, * BUFF2, BO, KBIND2, IRET) IF (IRET.NE.0) GO TO 990 C Subtract model from vis data. C Loop: Read vis. record. 200 CONTINUE ERRTXT = 'READING VIS FILE' CALL UVDISK ('READ', LUNI, FINDI, BUFF3, INIO2, IBIND, IRET) INIO = INIO2 IF (IRET.NE.0) GO TO 990 C Exit if no more data IF (INIO.LE.0) GO TO 300 C Uncompress Vis and put in AP CALL BUFPUT (APCORE, INIO, NRPARM, LREC, LLREC, IDATA, * BUFF3(IBIND), IBUFF(1)) C Loop over frequency groups C subtraction: normal LFR = 0 NFR = 0 IAPCC = IAPCC0 IF (.NOT.DOPTSP) THEN DO 220 I = I1,I2 IF (MODCHN(I).NE.LFR) THEN NFR = NFR + 1 LFR = MODCHN(I) IF (NFR.GT.1) IAPCC = IAPCC + MMCMP(NFR-1) * LENMOD MMCOMP = MMCMP(NFR) LCHAN = MODNCH(I) VIS = IDATA + LRPARM + (LFR-1)*JNCF + (VOFF/INCS)*JNCS C Fill frequency table C Double array to force accuracy FFRAC = (FREQG(LFR) / FREQ) - 1.0D0 FBUFF(1) = FFRAC IF (LCHAN.GT.1) THEN DO 210 J = 2,LCHAN INDEX = LFR + J - 1 FFRAC = (FREQG(INDEX) / FREQ) - 1.0D0 FBUFF(J) = FFRAC 210 CONTINUE END IF CALL QWD CALL QPUT (APCORE, FBUFF, 2, MCHAN, 2) CALL QWAIT C Do the arithmetic: subtract C Point IF (MODTYP.EQ.0) CALL QXXPTS (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, * NNCOR, SFLAG) C Gaussian IF (MODTYP.EQ.1) CALL QGASUB (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, * NNCOR, SFLAG) C Sphere IF (MODTYP.EQ.3) CALL QSPSUB (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, * NNCOR, SFLAG) END IF 220 CONTINUE ELSE LCHAN = 1 MMCOMP = 1 DO 230 LFR = CHANEL,CHANEL+NCHAN-1 X = FREQG(LFR) / FREQ X = LOG10 (X) Y = X*PTSPIX(1) + X*X*PTSPIX(2) + X*X*X*PTSPIX(3) + * X*X*X*X*PTSPIX(4) BTEMP = PTFL0 * (10.0D0 ** Y) CALL QPUT (APCORE, BTEMP, IAPCC, 1, 2) CALL QWD FFRAC = (FREQG(LFR) / FREQ) - 1.0D0 FBUFF(1) = FFRAC CALL QPUT (FBUFF, 2, 1, 2) CALL QWAIT VIS = IDATA + LRPARM + (LFR-1)*JNCF + (VOFF/INCS)*JNCS C Do the arithmetic: subtract C Point IF (MODTYP.EQ.0) CALL QXXPTS (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, NNCOR, * SFLAG) C Gaussian IF (MODTYP.EQ.1) CALL QGASUB (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, NNCOR, * SFLAG) C Sphere IF (MODTYP.EQ.3) CALL QSPSUB (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, MCHAN, NNCOR, * SFLAG) 230 CONTINUE END IF CALL QWR C Get UVs from AP (maybe pack UV) CALL BUFGET (APCORE, INIO, NRPARM, LREC, LLREC, IDATA, * BUFF1(KBIND1), IBUFF(1)) C Write vis record. NIOUT = INIO ERRTXT = 'WRITING VIS FILE' CALL UVDISK ('WRIT', LUN1, FIND1, BUFF1, NIOUT, KBIND1, IRET) IF (IRET.NE.0) GO TO 990 C Uncompress Vis and put in AP CALL BUFPUT (APCORE, INIO, NRPARM, LREC, LLREC, IDATA, * BUFF3(IBIND), IBUFF(1)) C Loop over frequency groups C division LFR = 0 NFR = 0 IAPCC = IAPCC0 IF (.NOT.DOPTSP) THEN DO 260 I = I1,I2 IF (MODCHN(I).NE.LFR) THEN NFR = NFR + 1 LFR = MODCHN(I) IF (NFR.GT.1) IAPCC = IAPCC + MMCMP(NFR-1) * LENMOD MMCOMP = MMCMP(NFR) LCHAN = MODNCH(I) VIS = IDATA + LRPARM + (LFR-1)*JNCF + (VOFF/INCS)*JNCS C Fill frequency table C Double array to force accuracy FFRAC = (FREQG(LFR) / FREQ) - 1.0D0 FBUFF(1) = FFRAC IF (LCHAN.GT.1) THEN DO 250 J = 2,LCHAN INDEX = LFR + J - 1 FFRAC = (FREQG(INDEX) / FREQ) - 1.0D0 FBUFF(J) = FFRAC 250 CONTINUE END IF CALL QWD CALL QPUT (APCORE, FBUFF, 2, MCHAN, 2) CALL QWAIT C Do the arithmetic: divide C Point IF (MODTYP.EQ.0) CALL QPTDIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) C Gaussian IF (MODTYP.EQ.1) CALL QGADIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) C Sphere IF (MODTYP.EQ.3) CALL QSPDIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) END IF 260 CONTINUE ELSE LCHAN = 1 MMCOMP = 1 DO 270 LFR = CHANEL,CHANEL+NCHAN-1 X = FREQG(LFR) / FREQ X = LOG10 (X) Y = X*PTSPIX(1) + X*X*PTSPIX(2) + X*X*X*PTSPIX(3) + * X*X*X*X*PTSPIX(4) BTEMP = PTFL0 * (10.0D0 ** Y) CALL QPUT (APCORE, BTEMP, IAPCC, 1, 2) CALL QWD FFRAC = (FREQG(LFR) / FREQ) - 1.0D0 FBUFF(1) = FFRAC CALL QPUT (APCORE, FBUFF, 2, 1, 2) CALL QWAIT VIS = IDATA + LRPARM + (LFR-1)*JNCF + (VOFF/INCS)*JNCS C Do the arithmetic: subtract C Point IF (MODTYP.EQ.0) CALL QPTDIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) C Gaussian IF (MODTYP.EQ.1) CALL QGADIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) C Sphere IF (MODTYP.EQ.3) CALL QSPDIV (APCORE, IAPCC, UV, VIS, * LLREC, JNCF, JNCS, MMCOMP, INIO, LCHAN, NNCOR) 270 CONTINUE END IF CALL QWR C Get UVs from AP (maybe pack UV) CALL BUFGET (APCORE, INIO, NRPARM, LREC, LLREC, IDATA, * BUFF2(KBIND2), IBUFF(1)) C Write vis record. NIOUT = INIO ERRTXT = 'WRITING VIS FILE' CALL UVDISK ('WRIT', LUN2, FIND2, BUFF2, NIOUT, KBIND2, IRET) IF (IRET.NE.0) GO TO 990 C Check if time for % done C message. PCLST = PCLST + NIOUT IPCDNE = (100.0 / NTIMES) * (PCLST / PCTOT) + 0.5 IPCDNE = IPCDNE - MOD (IPCDNE, 10) C Write % done message. IF ((IPCDNE.GT.IPCLST) .AND. (DOMSG)) THEN WRITE (MSGTXT,1240) IPCDNE IF (IPCDNE.LE.100) CALL MSGWRT (2) IPCLST = IPCDNE END IF C Check if time for AP roller CALL QROLL (APCORE, IAPCT, IBUFF, JBUFSZ, IRET) IF (IRET.NE.0) GO TO 999 GO TO 200 C Finish write 300 CALL QRLSE NIOUT = 0 ERRTXT = 'FLUSHING VIS FILE' CALL UVDISK ('FLSH', LUN1, FIND1, BUFF1, NIOUT, KBIND1, IRET) IF (IRET.NE.0) GO TO 990 NIOUT = 0 CALL UVDISK ('FLSH', LUN2, FIND2, BUFF2, NIOUT, KBIND2, IRET) IF (IRET.NE.0) GO TO 990 C Update no. comp. left. XNCOMP = XNCOMP - MXCMP C Close files CALL ZCLOSE (LUNI, FINDI, IRET) CALL ZCLOSE (LUN1, FIND1, IRET) CALL ZCLOSE (LUN2, FIND2, IRET) IF (WESET) MODMAX = 0 IRET = 0 GO TO 999 C Error 990 WRITE(MSGTXT,1990,ERR=999) IRET, ERRTXT 995 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVADFT: MODEL TYPES =',I2,I2,' INCOMPATIBLE') 1002 FORMAT ('EVADFT: ILLEGAL MODEL TYPE =',I3) 1006 FORMAT ('EVADFT: fields',I5,' -',I5,' chns',I5,' -',I5,' in',I5, * ' CC models') 1070 FORMAT ('EVADFT WRONG MODTYP AT FIELD, COMP',I5,I9) 1090 FORMAT ('EVADFT: ERROR',I5,' READING CLEAN COMPS REC',I5) 1240 FORMAT ('Model computation is ',I5,' percent complete') 1990 FORMAT ('EVADFT: ERROR',I5,' ',A) END SUBROUTINE EVAUVS (IRET) C----------------------------------------------------------------------- C EVAUVS reads through the UV data sets created by the model division C and subtraction and does various statistics on them. C Output: C IRET I Error code C----------------------------------------------------------------------- INTEGER IRET C INCLUDE 'EVAUV.INC' INTEGER K, VOL, CNO, VO, BO, LUN, FIND, LENBU, INIO, BIND, IREC, * IPTR, NCOPY, I, J, NP, HLIM, LVIS, IR, IRM, L, NOUT(2), NNO, * NNT, NTOT(2), LMAX CHARACTER WHICH*1, PHNAME*48, UTYPE*2 LOGICAL T, F REAL VIS(3,MAXCIF), R, AVREAL(2), AVIMAG(2), AVAMP(2), * AVRRMS(2), AVIRMS(2), AVARMS(2), SUBTR, RMM(10), RM, AOUT(2), * AAO, RMB, RM3, RM3S(2) DOUBLE PRECISION AS, ASS, RS, RSS, IS, ISS, WTS, HST(UVHIST), * HSTS(UVHIST), HSTW(UVHIST), ASW, RMS, AV, AA, AAA, WT, A INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DCAT.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DSCD.INC' INCLUDE 'INCS:PSTD.INC' EQUIVALENCE (VIS, BUFF3) DATA BO, VO, LUN /1, 0, 24/ DATA T, F /.TRUE., .FALSE./ DATA RMM / 5.0, 4.5, 4.0, 3.5, 6*3.0/ C----------------------------------------------------------------------- RM3S(1) = APARM(1) RM3S(2) = APARM(2) LMAX = 7 RSCALE = MIN (CELLSG(1), CELLSG(2)) IF (RSCALE.EQ.0.0) RSCALE = MIN (ABS(XCELL(1)), ABS(XCELL(2))) / * 3600. IF (RSCALE.EQ.0.0) RSCALE = 0.1 / 3600.0 RSCALE = 1.0 / (DG2RAD * RSCALE) / (UVHIST - 2) IF (SC1FIL.GT.0) THEN VOL = SCRVOL(SC1FIL) CNO = SCRCNO(SC1FIL) UTYPE = 'SC' ELSE VOL = FVOL(-SC1FIL) CNO = FCNO(-SC1FIL) UTYPE = 'UV' END IF WHICH = '1' SUBTR = 0.0 CALL COPY (256, UVBLK, CATBLK) CALL UVPGET (IRET) COMPDT = CATBLK(KINAX).EQ.1 NCOPY = LREC - NRPARM IF (COMPDT) THEN LVIS = NCOPY ELSE LVIS = NCOPY / 3 END IF NP = MIN (NCOR, 2) HLIM = UVHIST C sum histograms etc IRM = 0 DO 100 K = 1,2 REMXN(1,K) = 1.E10 REMXN(2,K) = -1.E10 IMMXN(1,K) = 1.E10 IMMXN(2,K) = -1.E10 C Open input file. CALL ZPHFIL (UTYPE, VOL, CNO, 1, PHNAME, IRET) CALL ZOPEN (LUN, FIND, VOL, PHNAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'OPEN SCRATCH FILE', WHICH GO TO 990 END IF C Setup for init. RMS = 1.D5 AV = 0.0D0 DO 60 L = 1,LMAX LENBU = 0 CALL UVINIT ('READ', LUN, FIND, NVIS, VO, LREC, LENBU, * JBUFSZ, BUFF1, BO, BIND, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'INIT SCRATCH IO', WHICH GO TO 990 END IF C zero counters RM = RMM(L) * RMS RMB = 6.0 * RMS NNT = 0 WTS = 0.0D0 AA = 0.0D0 AAA = 0.0D0 IF (L.EQ.LMAX) THEN RM3 = RM3S(K) * RMS CALL DFILL (HLIM, 0.0D0, HST) CALL DFILL (HLIM, 0.0D0, HSTS) CALL DFILL (HLIM, 0.0D0, HSTW) RS = 0.0D0 RSS = 0.0D0 IS = 0.0D0 ISS = 0.0D0 NNO = 0 AAO = 0 END IF C read loop 20 CALL UVDISK ('READ', LUN, FIND, BUFF1, INIO, BIND, IRET) IF (IRET.GT.0) THEN WRITE (MSGTXT,1000) IRET, 'READ SCRATCH DATA', WHICH GO TO 990 ELSE IF (INIO.GT.0) THEN IPTR = BIND DO 50 IREC = 1,INIO R = SQRT (BUFF1(IPTR+ILOCU)**2 + BUFF1(IPTR+ILOCV)**2) IR = (R / RSCALE) + 2.0 IR = MAX (1, MIN (HLIM, IR)) IF (COMPDT) THEN CALL ZUVXPN (NCOPY, BUFF1(IPTR+NRPARM), * BUFF1(IPTR+WTLOC), VIS) ELSE CALL RCOPY (NCOPY, BUFF1(IPTR+NRPARM), VIS) END IF DO 40 J = 1,NP DO 30 I = J,LVIS,NCOR WT = VIS(3,I) IF (WT.GT.0.0) THEN NNT = NNT + 1 VIS(1,I) = VIS(1,I) - SUBTR A = SQRT (VIS(1,I)**2 + VIS(2,I)**2) C accumulate for averages IF (ABS(A-AV).LE.RM) THEN AA = AA + WT * A AAA = AAA + WT * A * A WTS = WTS + WT END IF C get range values for Re/Im plots IF (L.EQ.LMAX) THEN IF (ABS(A-AV).LE.RM3) THEN REMXN(1,K) = MIN (REMXN(1,K), VIS(1,I)) REMXN(2,K) = MAX (REMXN(2,K), VIS(1,I)) IMMXN(1,K) = MIN (IMMXN(1,K), VIS(2,I)) IMMXN(2,K) = MAX (IMMXN(2,K), VIS(2,I)) END IF C accumulate histogram IF (ABS(A-AV).LE.RM) THEN HST(IR) = HST(IR) + WT * A HSTS(IR) = HSTS(IR) + WT * A * A HSTW(IR) = HSTW(IR) + WT RS = RS + WT * VIS(1,I) RSS = RSS + WT * VIS(1,I)**2 IS = IS + WT * VIS(2,I) ISS = ISS + WT * VIS(2,I)**2 C count really bad ones ELSE IF (ABS(A-AV).GT.RMB) THEN NNO = NNO + 1 AAO = AAO + ABS (A-AV) END IF END IF END IF 30 CONTINUE 40 CONTINUE IPTR = IPTR + LREC 50 CONTINUE GO TO 20 END IF IF (WTS.GT.0) THEN AV = AA / WTS AAA = AAA / WTS - AV * AV RMS = SQRT (MAX (0.0D0, AAA)) END IF 60 CONTINUE C close CALL ZCLOSE (LUN, FIND, IRET) C average counters AS = 0.0D0 ASS = 0.0D0 ASW = 0.0D0 DO 90 IR = 1,HLIM IF (HSTW(IR).GT.0.0) THEN AS = AS + HST(IR) ASS = ASS + HSTS(IR) ASW = ASW + HSTW(IR) HST(IR) = HST(IR) / HSTW(IR) HSTS(IR) = HSTS(IR) / HSTW(IR) - HST(IR)**2 HISTR(IR,K) = SQRT (MAX (0.0D0, HSTS(IR))) HISTA(IR,K) = HST(IR) IRM = MAX (IRM, IR) ELSE HISTA(IR,K) = 0.0 HISTR(IR,K) = 0.0 END IF 90 CONTINUE IF (ASW.GT.0.0) THEN AS = AS / ASW ASS = ASS / ASW - AS * AS ASS = SQRT (MAX (0.0D0, ASS)) END IF IF (WTS.GT.0.0) THEN RS = RS / WTS RSS = RSS / WTS - RS * RS RSS = SQRT (MAX (0.0D0, RSS)) IS = IS / WTS ISS = ISS / WTS - IS * IS ISS = SQRT (MAX (0.0D0, ISS)) END IF AVREAL(K) = RS AVRRMS(K) = RSS AVIMAG(K) = IS AVIRMS(K) = ISS AVAMP(K) = AS AVARMS(K) = ASS NOUT(K) = NNO AOUT(K) = AAO / MAX (1, NNO) NTOT(K) = NNT C point to second data set IF (SC2FIL.GT.0) THEN VOL = SCRVOL(SC2FIL) CNO = SCRCNO(SC2FIL) UTYPE = 'SC' ELSE VOL = FVOL(-SC2FIL) CNO = FCNO(-SC2FIL) UTYPE = 'UV' END IF WHICH = '2' SUBTR = 1.0 100 CONTINUE C report results WRITE (MSGTXT,1100) CALL MSGWRT (5) WRITE (MSGTXT,1101) 'subtract', AVREAL(1), AVRRMS(1), AVIMAG(1), * AVIRMS(1), AVAMP(1), AVARMS(1) CALL MSGWRT (5) WRITE (MSGTXT,1101) 'divide-1', AVREAL(2), AVRRMS(2), AVIMAG(2), * AVIRMS(2), AVAMP(2), AVARMS(2) CALL MSGWRT (5) WRITE (MSGTXT,1110) CALL MSGWRT (5) WRITE (MSGTXT,1111) 'subtract', NOUT(1), NTOT(1), AOUT(1) CALL MSGWRT (5) WRITE (MSGTXT,1111) 'divide-1', NOUT(2), NTOT(2), AOUT(2) CALL MSGWRT (5) RPARM(1) = AVREAL(1) RPARM(2) = AVRRMS(1) RPARM(3) = AVIMAG(1) RPARM(4) = AVIRMS(1) RPARM(5) = AVAMP(1) RPARM(6) = AVARMS(1) RPARM(7) = FLOAT (NOUT(1)) / FLOAT (NTOT(1)) RPARM(8) = AOUT(1) RPARM(9) = AVREAL(2) RPARM(10) = AVRRMS(2) RPARM(11) = AVIMAG(2) RPARM(12) = AVIRMS(2) RPARM(13) = AVAMP(2) RPARM(14) = AVARMS(2) RPARM(15) = FLOAT (NOUT(2)) / FLOAT (NTOT(2)) RPARM(16) = AOUT(2) GO TO 999 C 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVAUVS: ERROR',I4,' ON ',A,1X,A) 1100 FORMAT ('method',4X,'real part',10X,'imaginary part',5X, * 'amplitude') 1101 FORMAT (A8,3(F8.4,' +-',F7.3,1X)) 1110 FORMAT ('method',3X,'# bad samples',2X,'total samples',2X, * 'avg bad amp') 1111 FORMAT (A8,I14,I15,F13.4) END SUBROUTINE EVAIMG (NGAUSS, MAHIST, PXRANG, MSCALE, MUNITS, HISTM, * MFIT, RPARM, BUFF, IRET) C----------------------------------------------------------------------- C Uses the CCCNO, CCDISK to report statistics on the model images C Inputs: C MAHIST I Size of histogram C NGAUSS I Number of Gaussians in MFIELD C PXRANG R(2) User specified range for 1st histogram C Outputs: C MSCALE R(2,2) Max, min of total, rms-sized histograms C MUNITS C*8 Image units C HISTM I(*,2) Total range histogram, raange near avg C MFIT R(2) Robust average and rms C BUFF R(*) scratch buffer for image IO C RPARM R(30) returned adverb array C IRET I Error C----------------------------------------------------------------------- INTEGER NGAUSS, MAHIST, HISTM(MAHIST,2), IRET REAL PXRANG(2), MSCALE(2,2), MFIT(2), BUFF(*), RPARM(30) CHARACTER MUNITS*8 C INTEGER CATCLN(256), IFIELD, NFIELD, SCRTCH(256), NX, NY, * NWORDS, MLUN, MIND, NSUM LONGINT MPOFF REAL CATCR(256), MP(2) DOUBLE PRECISION AV, RMS, RSUM, RSUMS, ASUM, ASUMS CHARACTER STAT*4, PHNAME*48 LOGICAL FIRST INCLUDE 'INCS:PUVD.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DHDR.INC' EQUIVALENCE (CATCLN, CATCR) DATA MLUN /29/ C----------------------------------------------------------------------- NFIELD = MFIELD / MAX (1, NGAUSS) RSUM = 0.0 RSUMS = 0.0 ASUM = 0.0 ASUMS = 0.0 NSUM = 0 CELLSG(1) = 1.E6 CELLSG(2) = 1.E6 FIRST = .TRUE. CALL FILL (2*MAHIST, 0, HISTM) DO 100 IFIELD = 1,NFIELD STAT = 'REST' CALL CATIO ('READ', CCDISK(IFIELD), CCCNO(IFIELD), CATCLN, * STAT, SCRTCH, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IFIELD, IRET, 'READ CATBLK' GO TO 980 END IF IF (CATCR(KRCIC).NE.0.0) CELLSG(1) = MIN (CELLSG(1), * ABS(CATCR(KRCIC))) IF (CATCR(KRCIC+1).NE.0.0) CELLSG(2) = MIN (CELLSG(2), * ABS(CATCR(KRCIC+1))) NX = CATCLN(KINAX) NY = CATCLN(KINAX+1) IF (IFIELD.EQ.1) CALL H2CHR (8, 1, CATCR(KHBUN), MUNITS) C get memory NWORDS = (NX * NY - 1) / 1024 + 3 CALL ZMEMRY ('GET ', 'EVAIMG', NWORDS, MP, MPOFF, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IFIELD, IRET, 'GETTING MEMORY' GO TO 980 END IF C opn map file CALL ZPHFIL ('MA', CCDISK(IFIELD), CCCNO(IFIELD), 1, PHNAME, * IRET) CALL ZOPEN (MLUN, MIND, CCDISK(IFIELD), PHNAME, .TRUE., .TRUE., * .TRUE., IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IFIELD, IRET, 'OPENING MA FILE' GO TO 980 END IF CALL EVARMS (MLUN, MIND, NX, NY, MP(1+MPOFF), AV, RMS, BUFF, * IRET) IF (IRET.GT.0) GO TO 999 IF (IRET.EQ.0) THEN ASUM = ASUM + AV ASUMS = ASUMS + AV * AV RSUM = RSUM + RMS RSUMS = RSUMS + RMS * RMS NSUM = NSUM + 1 IF ((FIRST) .AND. (RMS.GT.0.0)) THEN IF (PXRANG(2).GT.PXRANG(1)) THEN MSCALE(1,1) = PXRANG(1) MSCALE(2,1) = PXRANG(2) ELSE MSCALE(1,1) = CATCR(KRDMN) MSCALE(2,1) = CATCR(KRDMX) END IF MSCALE(1,2) = AV - 5.*RMS MSCALE(2,2) = AV + 5.*RMS FIRST = .FALSE. END IF CALL EVAHST (NX, NY, MP(1+MPOFF), MAHIST, MSCALE, HISTM) END IF CALL ZCLOSE (MLUN, MIND, IRET) CALL ZMEMRY ('FREE', 'EVAIMG', NWORDS, MP, MPOFF, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IFIELD, IRET, 'GETTING MEMORY' GO TO 980 END IF 100 CONTINUE C average the averages/rmses IF (NSUM.GT.0) THEN ASUM = ASUM / NSUM ASUMS = ASUMS / NSUM - ASUM * ASUM ASUMS = SQRT (MAX (0.0D0, ASUMS)) RSUM = RSUM / NSUM RSUMS = RSUMS / NSUM - RSUM * RSUM RSUMS = SQRT (MAX (0.0D0, RSUMS)) IF (NFIELD.GT.1) THEN WRITE (MSGTXT,1100) ASUM, ASUMS, NFIELD CALL MSGWRT (5) WRITE (MSGTXT,1101) RSUM, RSUMS, NFIELD CALL MSGWRT (5) RPARM(18) = ASUMS RPARM(20) = RSUMS ELSE WRITE (MSGTXT,1110) ASUM, RSUM CALL MSGWRT (5) END IF MFIT(1) = ASUM MFIT(2) = RSUM RPARM(17) = ASUM RPARM(19) = RSUM ELSE MSGTXT = 'EVAIMG DID NOT FIND ANY IMAGE RMS TO REPORT' CALL MSGWRT (8) MFIT(1) = 0.0 MFIT(2) = 0.0 END IF IRET = 0 GO TO 999 C 980 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVAIMG FIELD',I5,' ERROR',I3,' ON ',A) 1100 FORMAT ('Average mean',1PE12.4,' +-',1PE11.4,' over',I5,' images') 1101 FORMAT ('Average rms ',1PE12.4,' +-',1PE11.4,' over',I5,' images') 1110 FORMAT ('Image mean',1PE12.4,' rms',1PE11.4) END SUBROUTINE EVARMS (MLUN, MIND, NX, NY, MP, AV, RMS, BUFF, IRET) C----------------------------------------------------------------------- C EVARMS does the robust rms computation on an image C Inputs: C MLUN I LUN of open image file C MIND I FTAB pointer for the file C NX I Number columns in image C NY I Number rows in image C Outputs: C MP R(*) Image-sized work area C AV R Robust mean of image C RMS R Robust average of image C BUFF R(*) IO buffer C IRET I Error code: -1 => no points found C----------------------------------------------------------------------- INTEGER MLUN, MIND, NX, NY, IRET REAL MP(NX,NY), BUFF(*) DOUBLE PRECISION AV, RMS C INTEGER IX, IY, WIN(4), NBY, IBLKOF, L, BPOS, NSUM REAL RM(7) DOUBLE PRECISION SUM, SUMS, RMM, V INCLUDE 'INCS:ZPBUFSZ.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DMSG.INC' DATA RM /6., 5., 4., 4*3./ C----------------------------------------------------------------------- C init IO NBY = 2 * UVBFSS WIN(1) = 1 WIN(2) = 1 WIN(3) = NX WIN(4) = NY IBLKOF = 1 CALL MINIT ('READ', MLUN, MIND, NX, NY, WIN, BUFF, NBY, IBLKOF, * IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'INIT IMAGE IO' GO TO 990 END IF C read in image DO 10 IY = 1,NY CALL MDISK ('READ', MLUN, MIND, BUFF, BPOS, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'READ IMAGE ROWS' GO TO 990 END IF CALL RCOPY (NX, BUFF(BPOS), MP(1,IY)) 10 CONTINUE C robust average RMS = 1.D6 AV = 0 DO 50 L = 1,7 SUM = 0.0D0 SUMS = 0.0D0 NSUM = 0 RMM = RM(L) * RMS DO 30 IY = 1,NY DO 20 IX = 1,NX IF (MP(IX,IY).NE.FBLANK) THEN V = MP(IX,IY) IF (ABS(V-AV).LE.RMM) THEN SUM = SUM + V SUMS = SUMS + V * V NSUM = NSUM + 1 END IF END IF 20 CONTINUE 30 CONTINUE IF (NSUM.GT.0.0D0) THEN SUM = SUM / NSUM SUMS = SUMS / NSUM - SUM * SUM SUMS = SQRT (MAX (0.0D0, SUMS)) AV = SUM RMS = SUMS ELSE IRET = -1 GO TO 999 END IF 50 CONTINUE GO TO 999 C 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVARMS: ERROR',I4,' ON ',A) END SUBROUTINE EVAHST (NX, NY, MP, MAHIST, MSCALE, HISTM) C----------------------------------------------------------------------- C EVAHST takes the image in ram and adds it to the histograms C Inputs: C NX I Number columns in image C NY I Number rows in image C MP R(*,*) Image data C MAHIST I Size of histogram + edge counters C MSCALE R(2,2) min/max for 2 histograms C In/out: C HISTM I(*,2) histogram C----------------------------------------------------------------------- INTEGER NX, NY, MAHIST, HISTM(MAHIST,2) REAL MP(NX,NY), MSCALE(2,2) C INTEGER IX, IY, I REAL DELTA1, DELTA2 INCLUDE 'INCS:DDCH.INC' C----------------------------------------------------------------------- IF (MSCALE(2,1).GT.MSCALE(1,1)) THEN DELTA1 = (MSCALE(2,1) - MSCALE(1,1)) / (MAHIST-2) ELSE DELTA1 = 1.0 END IF IF (MSCALE(2,2).GT.MSCALE(1,2)) THEN DELTA2 = (MSCALE(2,2) - MSCALE(1,2)) / (MAHIST-2) ELSE DELTA2 = 1.0 END IF DO 20 IY = 1,NY DO 10 IX = 1,NX IF (MP(IX,IY).NE.FBLANK) THEN I = (MP(IX,IY) - MSCALE(1,1)) / DELTA1 IF (MP(IX,IY).EQ.MSCALE(2,1)) I = MIN (I, MAHIST-3) I = MIN (MAHIST, MAX (1, I+2)) HISTM(I,1) = HISTM(I,1) + 1 I = (MP(IX,IY) - MSCALE(1,2)) / DELTA2 IF (MP(IX,IY).EQ.MSCALE(2,2)) I = MIN (I, MAHIST-3) I = MIN (MAHIST, MAX (1, I+2)) HISTM(I,2) = HISTM(I,2) + 1 END IF 10 CONTINUE 20 CONTINUE C 999 RETURN END SUBROUTINE EVAPLT (NP, MAP, SC1MAP, SC2MAP, IRET) C----------------------------------------------------------------------- C EVAPLT plots the histograms C Input: C NP I Size of work space C Output: C MAP D(*) work space C IRET I Error code C----------------------------------------------------------------------- INTEGER NP, IRET REAL SC1MAP(NP,*), SC2MAP(NP,*) DOUBLE PRECISION MAP(NP,*) C INCLUDE 'EVAUV.INC' INTEGER K, IM, IVER, GRCHN, TVCHN, TVCORN(2), PLUN, PIND, LABEL, * IROUND, I, VOL, CNO, LUN, FIND, LENBU, VO, BO, BIND, INIO, * IPTR, IREC, NCOPY, LVIS, NPP, WTLOC, J, IX, IY, NSM, XPRM REAL YMIN, YMAX, XMIN, XMAX, MHIST(MAHIST), BLOG, TEMP, DX, * DY, VIS(3,MAXCIF), X, Y, WT, SUBTR LOGICAL DOTV, DOLOG, T, F, COMPDT DOUBLE PRECISION PARMS(3), NGOOD, NBAD CHARACTER PFILE*48, UUNITS*25, WHICH*1, PHNAME*48, UTYPE*2 INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DGPH.INC' INCLUDE 'INCS:DCAT.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DUVH.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DGDS.INC' EQUIVALENCE (VIS, BUFF3) DATA T, F /.TRUE., .FALSE./ DATA LUN /24/ DATA VO, BO /0, 1/ C----------------------------------------------------------------------- IF (DOHIST.LE.0.0) GO TO 999 NSM = APARM(4) + 0.1 DOTV = XDOTV.GT.0.0 CALL FNDEXT ('PL', UVBLK, IVER) TVCHN = 1 TVCORN(1) = 0 TVCORN(2) = 0 GRCHN = XGRCH + 0.5 LABEL = IROUND (XLABEL) IF (MOD(LABEL,100).EQ.0) LABEL = (LABEL/100)*100 + 3 DOLOG = .FALSE. PARMS(1) = 0.0D0 XPRM = NPARM + 10 IF (DOHIST.GE.1.5) THEN UUNITS = 'Baseline length (lambdas)' DO 20 K = 1,2 IM = 0 YMIN = 1.E10 YMAX = -1.E10 DO 10 I = 1,UVHIST IF (HISTA(I,K).GT.0.0) THEN IM = I YMIN = MIN (YMIN, HISTA(I,K)-HISTR(I,K)) YMAX = MAX (YMAX, HISTA(I,K)+HISTR(I,K)) END IF 10 CONTINUE IF (IM.LT.4) THEN WRITE (MSGTXT,1010) IM, K CALL MSGWRT (6) GO TO 20 END IF IF ((YMIN.GT.0.0) .AND. (YMIN.LT.0.33*YMAX)) YMIN = 0.0 IVER = IVER + 1 IF (.NOT.DOTV) THEN CALL MADDEX ('PL', DISKIN, OLDCNO, INCATB, BUFF2, .TRUE., * 'UPDT', IVER, IRET) IF (IRET.NE.0) THEN MSGTXT = 'ERROR ADDING PLOT TO HEADER' CALL MSGWRT (7) END IF END IF IM = IM + 1 HISTA(IM,K) = HISTA(UVHIST,K) XMIN = 0.0 XMAX = IM * RSCALE XTRPRM(1) = K XTRPRM(2) = XMIN XTRPRM(3) = XMAX XTRPRM(4) = YMIN XTRPRM(5) = YMAX CALL ZPHFIL ('PL', DISKIN, OLDCNO, IVER, PFILE, IRET) CALL GINIT (DISKIN, OLDCNO, PFILE, 0, 39, XPRM, XNAMEI, * DOTV, TVCHN, GRCHN, TVCORN, INCATB, BUFF1, PLUN, PIND, * IRET) IF (IRET.NE.0) GO TO 20 CALL HISTOG (1+K, IM, HISTA(1,K), HISTR(1,K), XMIN, XMAX, * YMIN, YMAX, DOLOG, LABEL, IVER, PARMS, UUNITS, IBUFF1, * IRET) IF (IRET.NE.0) THEN MSGTXT = 'ERROR. WILL TRY TO FINISH PARTIAL GRAPH.' CALL MSGWRT (7) END IF GPHPAG = .TRUE. CALL GFINIS (BUFF1, IRET) IF (IRET.LT.0) GO TO 999 C Successful plot file finished. IF (IRET.EQ.0) THEN IF (.NOT.DOTV) THEN CALL HIPLOT (DISKIN, OLDCNO, IVER, BUFF2, IRET) WRITE (MSGTXT,1015) IVER CALL MSGWRT (5) IRET = 0 END IF END IF 20 CONTINUE C plot image histograms IF (.NOT.DOPTMD) THEN DOLOG = FUNCTY.EQ.'LG' BLOG = ALOG10 (0.5) PARMS(1) = 1.0D0 PARMS(2) = MFIT(1) PARMS(3) = MFIT(2) DO 40 K = 1,2 YMIN = 1.E10 YMAX = -1.E10 MHIST(1) = HISTM(1,K) MHIST(MAHIST) = HISTM(MAHIST,K) DO 30 I = 2,MAHIST-1 IF (DOLOG) THEN IF (HISTM(I,K).GT.0) THEN MHIST(I) = ALOG10 (HISTM(I,K)+0.0) ELSE MHIST(I) = BLOG END IF ELSE MHIST(I) = HISTM(I,K) END IF YMIN = MIN (YMIN, MHIST(I)) YMAX = MAX (YMAX, MHIST(I)) 30 CONTINUE IF ((YMIN.GT.0.0) .AND. (YMIN.LT.0.33*YMAX)) YMIN = 0.0 IVER = IVER + 1 IF (.NOT.DOTV) THEN CALL MADDEX ('PL', DISKIN, OLDCNO, INCATB, BUFF2, * .TRUE., 'UPDT', IVER, IRET) IF (IRET.NE.0) THEN MSGTXT = 'ERROR ADDING PLOT TO HEADER' CALL MSGWRT (7) END IF END IF XMIN = MSCALE(1,K) XMAX = MSCALE(2,K) XTRPRM(1) = K + 2.0 XTRPRM(2) = XMIN XTRPRM(3) = XMAX XTRPRM(4) = YMIN XTRPRM(5) = YMAX CALL ZPHFIL ('PL', DISKIN, OLDCNO, IVER, PFILE, IRET) CALL GINIT (DISKIN, OLDCNO, PFILE, 0, 39, XPRM, XNAMEI, * DOTV, TVCHN, GRCHN, TVCORN, INCATB, BUFF1, PLUN, PIND, * IRET) IF (IRET.NE.0) GO TO 40 CALL HISTOG (1, MAHIST, MHIST, MHIST, XMIN, XMAX, YMIN, * YMAX, DOLOG, LABEL, IVER, PARMS, MUNITS, IBUFF1, IRET) IF (IRET.NE.0) THEN MSGTXT = 'ERROR. WILL TRY TO FINISH PARTIAL GRAPH.' CALL MSGWRT (7) END IF GPHPAG = .TRUE. CALL GFINIS (BUFF1, IRET) IF (IRET.LT.0) GO TO 999 C Successful plot file finished. IF (IRET.EQ.0) THEN IF (.NOT.DOTV) THEN CALL HIPLOT (DISKIN, OLDCNO, IVER, BUFF2, IRET) WRITE (MSGTXT,1015) IVER CALL MSGWRT (5) IRET = 0 END IF END IF 40 CONTINUE END IF END IF C plot Real vs Imag IF (SC1FIL.GT.0) THEN VOL = SCRVOL(SC1FIL) CNO = SCRCNO(SC1FIL) UTYPE = 'SC' ELSE VOL = FVOL(-SC1FIL) CNO = FCNO(-SC1FIL) UTYPE = 'UV' END IF WHICH = '1' CALL COPY (256, UVBLK, CATBLK) COMPDT = CATBLK(KINAX).EQ.1 CALL UVPGET (IRET) NCOPY = LREC - NRPARM IF (COMPDT) THEN LVIS = NCOPY ELSE LVIS = NCOPY / 3 END IF NPP = MIN (NCOR, 2) IF (COMPDT) THEN CALL AXEFND (8, 'WEIGHT ', CATBLK(KIPCN), CATH(KHPTP), * WTLOC, IRET) IF ((IRET.NE.0) .OR. (WTLOC.LT.0)) THEN IRET = 5 MSGTXT = 'CANNOT FIND WEIGHT AND SCALE FOR COMPRESSED DATA' CALL MSGWRT (8) GO TO 999 END IF IRET = 0 END IF SUBTR = 0.0 DO 200 K = 1,2 I = NP * NP CALL DFILL (I, 0.0D0, MAP) YMIN = IMMXN(1,K) YMAX = IMMXN(2,K) IF ((YMIN.GT.0.0) .AND. (YMIN.LT.0.33*YMAX)) YMIN = 0.0 TEMP = YMAX - YMIN IF (ABS(YMAX+YMIN).LT.0.05*TEMP) THEN TEMP = MAX (YMAX, -YMIN) YMAX = TEMP YMIN = -TEMP TEMP = YMAX - YMIN END IF YMAX = YMAX + 0.08 * TEMP YMIN = YMIN - 0.08 * TEMP XMIN = REMXN(1,K) XMAX = REMXN(2,K) IF ((XMIN.GT.0.0) .AND. (XMIN.LT.0.33*XMAX)) XMIN = 0.0 TEMP = XMAX - XMIN DX = TEMP / 400.0 * 1000.0 DY = TEMP / 400.0 * 1000.0 IF (ABS(XMAX+XMIN).LT.0.05*TEMP) THEN TEMP = MAX (XMAX, -XMIN) XMAX = TEMP XMIN = -TEMP TEMP = XMAX - XMIN END IF XMAX = XMAX + 0.08 * TEMP XMIN = XMIN - 0.08 * TEMP IF ((ABS(XMAX-YMAX).LT.0.05*TEMP) .AND. * (ABS(XMIN-YMIN).LT.0.05*TEMP) .AND. * (ABS(XMAX+XMIN).LT.0.05*TEMP)) THEN XMAX = MAX (XMAX, YMAX) XMAX = MAX (XMAX, -XMIN) XMAX = MAX (XMAX, -YMIN) YMAX = XMAX XMIN = -XMAX YMIN = -XMAX END IF IVER = IVER + 1 IF (.NOT.DOTV) THEN CALL MADDEX ('PL', DISKIN, OLDCNO, INCATB, BUFF2, .TRUE., * 'UPDT', IVER, IRET) IF (IRET.NE.0) THEN MSGTXT = 'ERROR ADDING PLOT TO HEADER' CALL MSGWRT (7) END IF END IF XTRPRM(1) = K + 4.0 XTRPRM(2) = XMIN XTRPRM(3) = XMAX XTRPRM(4) = YMIN XTRPRM(5) = YMAX CALL ZPHFIL ('PL', DISKIN, OLDCNO, IVER, PFILE, IRET) CALL GINIT (DISKIN, OLDCNO, PFILE, 0, 39, XPRM, XNAMEI, DOTV, * TVCHN, GRCHN, TVCORN, INCATB, BUFF1, PLUN, PIND, IRET) IF (IRET.NE.0) GO TO 200 CALL REIMPL (K, XMIN, XMAX, YMIN, YMAX, LABEL, IVER, BUFF1, * IRET) IF (IRET.NE.0) GO TO 190 C Open input file. CALL ZPHFIL (UTYPE, VOL, CNO, 1, PHNAME, IRET) CALL ZOPEN (LUN, FIND, VOL, PHNAME, T, F, T, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'OPEN SCRATCH FILE', WHICH GO TO 990 END IF LENBU = 0 CALL UVINIT ('READ', LUN, FIND, NVIS, VO, LREC, LENBU, * JBUFSZ, BUFF2, BO, BIND, IRET) IF (IRET.NE.0) THEN WRITE (MSGTXT,1000) IRET, 'INIT SCRATCH IO', WHICH GO TO 990 END IF NGOOD = 0.0D0 NBAD = 0.0D0 100 CALL UVDISK ('READ', LUN, FIND, BUFF2, INIO, BIND, IRET) IF (IRET.GT.0) THEN WRITE (MSGTXT,1000) IRET, 'READ SCRATCH DATA', WHICH GO TO 990 ELSE IF (INIO.GT.0) THEN IPTR = BIND DO 130 IREC = 1,INIO IF (COMPDT) THEN CALL ZUVXPN (NCOPY, BUFF2(IPTR+NRPARM), * BUFF2(IPTR+WTLOC), VIS) ELSE CALL RCOPY (NCOPY, BUFF2(IPTR+NRPARM), VIS) END IF DO 120 J = 1,NPP DO 110 I = J,LVIS,NCOR WT = VIS(3,I) IF (WT.GT.0.0) THEN X = VIS(1,I) - SUBTR Y = VIS(2,I) IF ((X.LE.XMAX) .AND. (X.GE.XMIN) .AND. * (Y.LE.YMAX) .AND. (Y.GE.YMIN)) THEN NGOOD = NGOOD + 1.0D0 IX = (X - XMIN) / (XMAX - XMIN) * (NP - 1.) + * 1.5 IY = (Y - YMIN) / (YMAX - YMIN) * (NP - 1.) + * 1.5 MAP(IX,IY) = MAP(IX,IY) + 1.0D0 ELSE NBAD = NBAD + 1 END IF END IF 110 CONTINUE 120 CONTINUE IPTR = IPTR + LREC 130 CONTINUE GO TO 100 END IF C draw contours 190 IF (NGOOD.GT.0.0D0) CALL IMGDRW (NP, MAP, NSM, SC1MAP, SC2MAP, * XMIN, XMAX, YMIN, YMAX, IBUFF1, IRET) WRITE (MSGTXT,1190) NGOOD, NBAD, WHICH CALL MSGWRT (4) WRITE (MSGTXT,1191) NP, NSM CALL MSGWRT (4) GPHPAG = K.EQ.1 CALL GFINIS (BUFF1, IRET) IF (IRET.LT.0) GO TO 999 C Successful plot file finished. IF (IRET.EQ.0) THEN IF (.NOT.DOTV) THEN CALL HIPLOT (DISKIN, OLDCNO, IVER, BUFF2, IRET) WRITE (MSGTXT,1015) IVER CALL MSGWRT (5) IRET = 0 END IF END IF IF (SC2FIL.GT.0) THEN VOL = SCRVOL(SC2FIL) CNO = SCRCNO(SC2FIL) UTYPE = 'SC' ELSE VOL = FVOL(-SC2FIL) CNO = FCNO(-SC2FIL) UTYPE = 'UV' END IF WHICH = '2' CALL ZCLOSE (LUN, FIND, IRET) SUBTR = 1.0 200 CONTINUE GO TO 999 C 990 CALL MSGWRT (8) C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('ERROR',I5,' ON ',A,' ',A) 1010 FORMAT ('ONLY',I2,' BOXES SAMPLED IN UV HISTOGRAM',I2) 1015 FORMAT ('Successful histogram plot file version',I7,' created') 1190 FORMAT ('Plotted',F11.0,' omitted',F10.0,' vis of RE/IM plot ',A) 1191 FORMAT ('Used image of',I5,' pixels on a side, smoothed by',I3, * ' pixels') END SUBROUTINE HISTOG (ITY, NBOXES, HIST, HRMS, XMIN, XMAX, YMIN, * YMAX, DOLOG, LABEL, IVER, PARMS, UNITS, BUFF1, IRET) C----------------------------------------------------------------------- C This routine will write commands to an open plot file for drawing C a histogram. C Inputs: C ITY I Plot type: 1 Flux on X axis, count on Y axis C 2 Radius on X axis, Avg flux +- rms on Y C 3 Radius on X axis, Avg gain - (1,0) +- rms Y C NBOXES I number of boxes for histogram. C HIST R(*) Histogram (underflow in (1) and overflow in C (NBOXES)) C XMIN R Min value in X C XMAX R Max value in X C YMIN R Min value in Y C XMAX R Max value in Y C DOLOG L T => use log(n) rather than linear scale C LABEL I Type of labeling C IVER I Plot file version number C PARMS D(3) Fit parameters (0) = 0 => none C UNITS C*8 Units to use C In/out: C BUFF1 I(*) I/O buffer for open, initialized pl file. C Output: C IRET I error code. 0=ok, 1=write error to plot file. C----------------------------------------------------------------------- INTEGER ITY, NBOXES, LABEL, IVER, BUFF1(*), IRET REAL XMIN, XMAX, YMIN, YMAX, HIST(*), HRMS(*) DOUBLE PRECISION PARMS(3) LOGICAL DOLOG CHARACTER UNITS*(*) C REAL BLC(7), CH(4), TRC(7), X, Y, FAC, XYRATO, LOCRAN(2), * EXF, BFACT, AX(5), AY(5), RANGE(2), W, SF, SW, DCX, DCY, XDIST, * DX, DY INTEGER IDEPTH(5), I, LTYPE, JM CHARACTER TXTMSG*80 INCLUDE 'INCS:DMSG.INC' C----------------------------------------------------------------------- C Set character offsets. LOCRAN(2) = YMAX LOCRAN(1) = YMIN CALL GTICNT (LABEL, LOCRAN, I) RANGE(1) = XMIN RANGE(2) = XMAX C number characters around CALL RFILL (4, 0.5, CH) LTYPE = MOD (ABS (LABEL), 100) IF (LTYPE.EQ.2) CH(1) = 2.5 IF (LTYPE.GT.2) CH(1) = I + 4.0 IF (LTYPE.GT.1) CH(2) = 2.0 IF (LTYPE.GT.2) CH(2) = CH(2) + 1.333 IF ((LTYPE.GT.1) .AND. (LTYPE.LT.7)) CH(2) = CH(2) + 3 * 1.333 IF (LTYPE.EQ.2) CH(3) = 2.5 IF (LTYPE.GT.1) CH(4) = CH(4) + 1.5 IF (LTYPE.GT.2) CH(4) = CH(4) + 1.333 IF ((LTYPE.GT.1) .AND. (LTYPE.LT.7)) CH(4) = CH(4) + 1.333 IF ((LABEL.GT.1) .AND. (LTYPE.LT.7)) CH(4) = CH(4) + 1.333 C Set BLC, TRC, XYRATO. CALL RFILL (5, 1.0, BLC(3)) CALL RFILL (5, 1.0, TRC(3)) CALL FILL (5, 1, IDEPTH) BLC(1) = -1.0 TRC(1) = NBOXES + 1.0 BLC(2) = YMIN - 0.05 * (YMAX - YMIN) TRC(2) = YMAX + 0.05 * (YMAX - YMIN) XYRATO = (TRC(1) - BLC(1)) / (TRC(2) - BLC(2)) C Kludge to keep XYRATO in bounds C to prevent overflow in GINITL. FAC = 1.0 IF (XYRATO.GT.3.0) THEN DO 10 I = 1,10000 IF (XYRATO.LT.2.0) GO TO 20 FAC = FAC / 2. XYRATO = XYRATO / 2. 10 CONTINUE ELSE IF (XYRATO.LT.0.333) THEN DO 15 I = 1,10000 IF (XYRATO.GT.0.50) GO TO 20 FAC = FAC * 2. XYRATO = XYRATO * 2. 15 CONTINUE END IF C 20 TRC(2) = TRC(2) / FAC BLC(2) = BLC(2) / FAC XYRATO = 1.0 / XYRATO C Initialize plot file line drw. CALL GINITL (BLC, TRC, XYRATO, CH, IDEPTH, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GLTYPE (1, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Draw borders. CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (TRC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (TRC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Labeling. I = NBOXES - 2 CALL HLABX (ITY, BLC, TRC, FAC, I, RANGE, DOLOG, IVER, HIST(1), * HIST(NBOXES), LABEL, UNITS, BUFF1, XDIST, IRET) IF (IRET.NE.0) GO TO 999 TXTMSG = 'End labeling, draw histogram' CALL GCOMNT (-1, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GLTYPE (2, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Type 1 plot IF (ITY.EQ.1) THEN C Position at first data point. X = 1 Y = BLC(2) CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Loop for rest of data points. SW = 0.0 JM = 0 DO 30 I = 2,NBOXES-1 IF (HIST(I).GT.SW) THEN SW = HIST(I) JM = I END IF Y = HIST(I) / FAC CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 X = I CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 Y = BLC(2) CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 30 CONTINUE C Gaussian fit IF (PARMS(1).GT.0.0) THEN TXTMSG = 'Draw gaussian fit' CALL GCOMNT (2, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GLTYPE (4, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 JM = JM - 2 IF (JM.GE.0) THEN IF (JM.LT.NBOXES/3) THEN WRITE (TXTMSG,1030) PARMS(2) * XDIST X = TRC(1) Y = TRC(2) CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = -29.0 DCY = -3.0 CALL GCHAR (12, 0, DCX, DCY, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 WRITE (TXTMSG,1031) PARMS(3) * XDIST X = TRC(1) Y = TRC(2) CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = -14.0 CALL GCHAR (11, 0, DCX, DCY, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 ELSE WRITE (TXTMSG,1030) PARMS(2) * XDIST X = BLC(1) Y = TRC(2) CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = 3.0 DCY = -3.0 CALL GCHAR (12, 0, DCX, DCY, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 WRITE (TXTMSG,1031) PARMS(3) * XDIST X = TRC(1) Y = TRC(2) CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = -14.0 DCY = -3.0 CALL GCHAR (11, 0, DCX, DCY, TXTMSG, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF END IF C solve for peak SF = 0.0 SW = 0.0 IF (DOLOG) BFACT = 10.0 ** YMAX DO 35 I = 2,NBOXES-1 X = RANGE(1) + (I * (RANGE(2)-RANGE(1))) / NBOXES X = (X - PARMS(2)) / PARMS(3) X = X * X / 2.0 IF (X.LE.39.) THEN W = EXP(-X) IF ((DOLOG) .AND. (HIST(I).GE.0.0)) THEN SF = SF + (10.0**HIST(I)) / BFACT SW = SW + W ELSE IF (.NOT.DOLOG) THEN SF = SF + HIST(I) / YMAX SW = SW + W END IF END IF 35 CONTINUE IF (SW.GT.0.0) THEN IF (DOLOG) THEN W = SF / SW * BFACT W = LOG10 (W) ELSE W = SF / SW * YMAX END IF ELSE W = YMAX END IF EXF = NBOXES / 1000.0 BFACT = W IF (DOLOG) BFACT = 10.0 ** W DO 40 I = 1,1001 X = RANGE(1) + (I-1.0) * (RANGE(2)-RANGE(1)) / 1000.0 X = (X - PARMS(2)) / PARMS(3) X = X * X / 2.0 IF (X.LE.69.) THEN Y = BFACT * EXP (-X) ELSE Y = 0. END IF IF (DOLOG) Y = ALOG10 (MAX (0.5, Y)) X = (I - 1) * EXF Y = Y / FAC IF (I.EQ.1) THEN CALL GPOS (X, Y, BUFF1, IRET) ELSE CALL GVEC (X, Y, BUFF1, IRET) END IF IF (IRET.NE.0) GO TO 999 40 CONTINUE END IF C type 2 plot ELSE DO 50 I = 2,NBOXES-1 X = I - 1 Y = (HIST(I) - HRMS(I)) / FAC CALL GPOS (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 Y = (HIST(I) + HRMS(I)) / FAC CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 X = I CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 Y = (HIST(I) - HRMS(I)) / FAC CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 X = I - 1 CALL GVEC (X, Y, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 50 CONTINUE CALL GLTYPE (4, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DX = 0.4 DY = 0.4 IF (XYRATO.GT.1.0) THEN DY = DY * XYRATO ELSE DX = DX / XYRATO END IF DO 60 I = 2,NBOXES-1 AX(1) = I - 0.5 AY(1) = HIST(I) / FAC AX(2) = AX(1) AX(3) = AX(1) AX(4) = AX(1) - DX AX(5) = AX(1) + DX AY(2) = AY(1) + DY AY(3) = AY(1) - DY AY(4) = AY(1) AY(5) = AY(1) CALL PNTPLT (3, AX, AY, BLC, TRC, .FALSE., .FALSE., BUFF1, * IRET) IF (IRET.NE.0) GO TO 999 60 CONTINUE END IF C 999 RETURN C----------------------------------------------------------------------- 1030 FORMAT ('MEAN',F8.4) 1031 FORMAT ('RMS',F8.4) END SUBROUTINE HLABX (ITY, BLC, TRC, FAC, NBOXES, RANGE, DOLOG, IVER, * UNDER, OVER, LABEL, UNITS, BUFF1, XDIST, IRET) C----------------------------------------------------------------------- C Write labeling for histogram. C Inputs: C ITY I Type: 1 image, 2 vis, 3 gain C BLC R(2) bottom left corner of plot. C TRC R(2) top right hand corner of plot. C FAC R FAC*XYRATO = real XYRATIO. C IVER I plot file version number C LABEL I labeling type C UNITS C*8 units C In/out: C BUFF1 I(256) I/O buffer for plot file. C Output: C XDIST R Scaling of X axis C IRET I error code returned from GVEC. C----------------------------------------------------------------------- REAL BLC(7), TRC(7), FAC, RANGE(2), UNDER, OVER, XDIST INTEGER ITY, NBOXES, IVER, LABEL, BUFF1(256), IRET LOGICAL DOLOG CHARACTER UNITS*(*) C CHARACTER PREFIX*5, TIME*8, DATE*12, CTEMP*8, CSTOK(12)*4, * NAMSTR*18, MSGBUF*80 LOGICAL PFLAG REAL XINTER(21), DCX, DCY, XNOINT, DIST, ODIST, XMAX, TICSCL, * YTICEL, YTICER, XVAL, YPOS, TICLEN, XINT, X, FREQ, DCXM INTEGER INOINT, INCHAR, I, IXO, M, ITRY, NXFR, NXST, NAX, INC, * IANGL, JSTOK, IT(3), ID(3), ICPNT, ITMP INCLUDE 'INCS:DCAT.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DHDR.INC' DATA TICSCL /70.0/ DATA CSTOK /'????','Beam','Ipol','Qpol','Upol','Vpol','Ppol', * 'Fpol','Pang','Spix','Optd',' '/ DATA XINTER /.1, .2, .5, 1., 2., 5., 10., 20., 50., 100., 200., * 500., 1000., 2000., 5000., 10000., 20000., 50000., 100000., * 200000., 500000./ C----------------------------------------------------------------------- IF (MOD(ABS(LABEL),100).LE.1) GO TO 999 C Tic positions. TICLEN = (TRC(1) - BLC(1)) / TICSCL YTICEL = BLC(1) + TICLEN YTICER = TRC(1) - TICLEN C Find interval value. DIST = FAC * (TRC(2) - BLC(2)) XINT = 8.0 DO 20 I = 1,21 XNOINT = AINT (DIST/XINTER(I)) IF (XNOINT.LE.XINT) GO TO 30 20 CONTINUE GO TO 110 C Interval and no of inter found. 30 XINT = XINTER(I) INOINT = XNOINT + 2.5 XVAL = AINT (FAC*BLC(2)/XINT) * XINT IF (XVAL.GT.FAC*BLC(2)) XVAL = XVAL - XINT IXO = I DCXM = -0.5 C Loop for all tics. DO 100 I = 1,INOINT XVAL = XVAL + XINT YPOS = XVAL / FAC IF (YPOS.GT.TRC(2)) GO TO 110 C TOP tic. CALL GPOS (TRC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (YTICER, YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Left hand tic. CALL GPOS (YTICEL, YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Write value. IF (MOD (ABS(LABEL),100).GT.2) THEN WRITE (MSGBUF,1030) XVAL CALL CHTRIM (MSGBUF, 14, MSGBUF, INCHAR) IF (IXO.GT.3) INCHAR = INCHAR - 2 DCX = - INCHAR - 1.0 DCY = -0.5 DCXM = MIN (DCXM, DCX) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF 100 CONTINUE C Number of pixels 110 DCX = DCXM - 2.0 YPOS = (TRC(2) + BLC(2)) / 2.0 CALL GPOS (BLC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (ITY.EQ.1) THEN IF (DOLOG) THEN MSGBUF = 'Log10 (number of pixels)' INCHAR = 24 ELSE MSGBUF = 'Number of pixels' INCHAR = 16 END IF ELSE IF (ITY.EQ.2) THEN MSGBUF = 'Flux residual (Jy)' INCHAR = 18 ELSE MSGBUF = 'Gain - (1,0)' INCHAR = 12 END IF DCY = INCHAR / 2.0 - 1.0 CALL GCHAR (INCHAR, 1, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Write bucket numbers on top IF (NBOXES.LE.8) THEN M = 1 ELSE IF (NBOXES.LE.16) THEN M = 2 ELSE IF (NBOXES.LE.40) THEN M = 5 ELSE IF (NBOXES.LE.80) THEN M = 10 ELSE IF (NBOXES.LE.160) THEN M = 20 ELSE IF (NBOXES.LE.400) THEN M = 50 ELSE IF (NBOXES.LE.800) THEN M = 100 ELSE IF (NBOXES.LE.1600) THEN M = 200 END IF TICLEN = (TRC(2) - BLC(2)) / TICSCL YTICEL = BLC(2) + TICLEN YTICER = TRC(2) - TICLEN DCY = 0.5 DO 150 I = 0,NBOXES,M X = I + 0.5 CALL GPOS (X, YTICER, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (X, TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (MOD (ABS(LABEL),100).GT.2) THEN WRITE (MSGBUF,1115) I CALL CHTRIM (MSGBUF, 4, MSGBUF, INCHAR) DCX = 0.5 - REAL(INCHAR) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF 150 CONTINUE C Label RHS bucket # X = (TRC(1) + BLC(1)) / 2.0 CALL GPOS (X, TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = -5.0 DCY = 0.5 IF (MOD(ABS(LABEL),100).GT.2) DCY = DCY + 1.333 MSGBUF = 'Box number' CALL GCHAR (10, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Range = IF (MOD(ABS(LABEL),100).LT.7) THEN CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = 0.0 DCY = -2.833 IF (MOD(ABS(LABEL),100).GT.2) DCY = DCY - 1.333 WRITE (MSGBUF,1151) RANGE(1), RANGE(2), UNITS CALL REFRMT (MSGBUF, '_', INCHAR) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Interval = CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 X = (RANGE(2) - RANGE(1)) / NBOXES WRITE (MSGBUF,1152) X, UNITS DCY = DCY - 1.333 CALL REFRMT (MSGBUF, '_', INCHAR) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) C Underflow = overflow = CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 I = UNDER + 0.1 M = OVER + 0.1 WRITE (MSGBUF,1154) I, M CALL REFRMT (MSGBUF, '_', INCHAR) DCY = DCY - 1.333 CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) END IF C Determine label range DIST = RANGE(2) - RANGE(1) ODIST = DIST CALL METSCL (LABEL, DIST, PREFIX, PFLAG) IF (PFLAG) GO TO 190 XDIST = DIST / ODIST ODIST = XDIST * RANGE(1) C Get interval DO 160 ITRY = 1,21 XNOINT = AINT (DIST/XINTER(ITRY)) IF (XNOINT.LE.8.0) GO TO 170 160 CONTINUE GO TO 190 C Bottom (value) tics 170 XINT = XINTER(ITRY) DCY = -1.5 XMAX = MAX (ABS(RANGE(2)), ABS(RANGE(1))) * XDIST INOINT = XNOINT + 2.5 XVAL = AINT (ODIST/XINT) * XINT IF (XVAL.GE.ODIST) XVAL = XVAL - XINT DO 175 I = 1,INOINT XVAL = XVAL + XINT X = ((XVAL-ODIST)/DIST) * NBOXES + 1.0 IF (X.GT.TRC(1)) GO TO 180 CALL GPOS (X, YTICEL, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (X, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (MOD(ABS(LABEL),100).GT.2) THEN WRITE (MSGBUF,1030) XVAL CALL CHTRIM (MSGBUF, 14, MSGBUF, INCHAR) IF (ITRY.GT.3) INCHAR = INCHAR - 2 DCX = 0.5 - INCHAR CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF 175 CONTINUE C Label with prefix 180 DCY = -1.5 IF (MOD(ABS(LABEL),100).GT.2) DCY = -2.833 X = (TRC(1) + BLC(1)) / 2.0 CALL GPOS (X, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 WRITE (MSGBUF,1175) PREFIX, UNITS CALL CHTRIM (MSGBUF, 14, MSGBUF, INCHAR) DCX = 0.5 - INCHAR / 2.0 CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (MOD(ABS(LABEL),100).GE.7) GO TO 999 C which axis is which? 190 NXFR = 0 NXST = 0 NAX = CATBLK(KIDIM) INC = 2 DO 200 I = 1,NAX ICPNT = KHCTP+(I-1)*INC CALL H2CHR (8, 1, CATH(ICPNT), CTEMP) IF (CTEMP(1:4).EQ.'FREQ') NXFR = I IF (CTEMP(1:4).EQ.'STOK') NXST = I 200 CONTINUE C Source name, stokes, freq. CALL GPOS (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = 0.0 DCY = 1.833 IF (MOD(ABS(LABEL),100).GT.2) DCY = DCY + 1.333 IANGL = 0 CALL H2CHR (8, 1, CATH(KHOBJ), CTEMP) FREQ = 0.0 JSTOK = 12 IF (NXFR.GT.2) FREQ = CATD(KDCRV+NXFR-1) + CATR(KRCIC+NXFR-1) * * (BLC(NXFR) - CATR(KRCRP+NXFR-1)) FREQ = FREQ / 1.E6 IF (NXST.GT.2) JSTOK = CATD(KDCRV+NXST-1) + CATR(KRCIC+NXST-1) * * (BLC(NXST) - CATR(KRCRP+NXST-1)) + 2.5 IF (NXFR.GT.2) WRITE (MSGBUF,1200) CTEMP, CSTOK(JSTOK), * FREQ IF (NXFR.LE.2) WRITE (MSGBUF,1200) CTEMP, CSTOK(JSTOK) CALL REFRMT (MSGBUF, '_', INCHAR) C image name INCHAR = INCHAR + 1 IF (INCHAR.GT.1) THEN MSGBUF(INCHAR:INCHAR+2) = ' _' INCHAR = INCHAR + 3 END IF CALL H2CHR (12, KHIMNO, CATH(KHIMN), NAMSTR(1:12)) CALL H2CHR (6, KHIMCO, CATH(KHIMC), NAMSTR(13:18)) CALL NAMEST (NAMSTR, CATBLK(KIIMS), MSGBUF(INCHAR:), ITMP) CALL REFRMT (MSGBUF, '_', INCHAR) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C time/date, version IF (LABEL.GT.0) THEN CALL GPOS (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL ZDATE (ID) CALL ZTIME (IT) CALL TIMDAT (IT, ID, TIME, DATE) WRITE (MSGBUF,1210) IVER, DATE, TIME CALL REFRMT (MSGBUF, '_', INCHAR) DCY = DCY + 1.333 CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) END IF C 999 RETURN C----------------------------------------------------------------------- 1030 FORMAT (F14.1) 1115 FORMAT (I4) 1151 FORMAT ('Range =',1PE12.4,' to',1PE12.4,1X,A) 1152 FORMAT ('Interval =',1PE12.4,1X,A) 1154 FORMAT ('Underflow =',I10,' _Overflow =',I10) 1175 FORMAT (A5,1X,A) 1200 FORMAT (A,' _',A4,'_ ',F10.3,' MHz') 1210 FORMAT ('Plot file version',I4,'__created ',A,A) END SUBROUTINE REIMPL (K, XMIN, XMAX, YMIN, YMAX, LABEL, IVER, BUFF1, * IRET) C----------------------------------------------------------------------- C Does plot scale and other initialization, labels and returns C Inputs: C K I 1 - subtraction, 2 - division C XMIN R X minimum C XMAX R X maximum C YMIN R Y minimum C YMAX R Y maximum C LABEL I Labeling type code C IVER I Plot file version number C In/Out: C BUFF1 R(*) Plot buffer C Output: C IRET I Error code C----------------------------------------------------------------------- INTEGER K, LABEL, IVER, IRET REAL XMIN, XMAX, YMIN, YMAX, BUFF1(*) C REAL BLC(7), CH(4), TRC(7), X, FAC, XYRATO, LOCRAN(2), FREQ, * RANGE(2), DIST, ODIST, DCX, DCY, XINTER(21), XNOINT, TICLEN, * YTICEL, YTICER, XINT, XVAL, YPOS, DCXM INTEGER IDEPTH(5), I, LTYPE, INOINT, IXO, INCHAR, NXFR, NXST, * NAX, INC, ICPNT, IANGL, JSTOK, IT(3), ID(3), ITRY, ITMP CHARACTER PREFIX*5, TIME*8, DATE*12, CTEMP*8, CSTOK(12)*4, * NAMSTR*18, MSGBUF*80 LOGICAL PFLAG INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DCAT.INC' INCLUDE 'INCS:DHDR.INC' DATA XINTER /.1, .2, .5, 1., 2., 5., 10., 20., 50., 100., 200., * 500., 1000., 2000., 5000., 10000., 20000., 50000., 100000., * 200000., 500000./ C----------------------------------------------------------------------- C Scaling LOCRAN(2) = YMAX LOCRAN(1) = YMIN CALL GTICNT (LABEL, LOCRAN, I) RANGE(1) = XMIN RANGE(2) = XMAX DIST = MAX (LOCRAN(2)-LOCRAN(1), RANGE(2)-RANGE(1)) ODIST = DIST CALL METSCL (LABEL, DIST, PREFIX, PFLAG) FAC = DIST / ODIST C number characters around CALL RFILL (4, 0.5, CH) LTYPE = MOD (ABS (LABEL), 100) IF (LTYPE.EQ.2) CH(1) = 2.5 IF (LTYPE.GT.2) CH(1) = I + 4.0 IF (LTYPE.GT.1) CH(2) = 2.0 IF (LTYPE.GT.2) CH(2) = CH(2) + 1.333 IF ((LTYPE.GT.1) .AND. (LTYPE.LT.7)) CH(2) = CH(2) + 1.333 IF (LTYPE.EQ.2) CH(3) = 2.5 IF (LTYPE.GT.1) CH(4) = CH(4) + 1.5 IF ((LTYPE.GT.1) .AND. (LTYPE.LT.7)) CH(4) = CH(4) + 1.333 IF ((LABEL.GT.1) .AND. (LTYPE.LT.7)) CH(4) = CH(4) + 1.333 C Set BLC, TRC, XYRATO. CALL RFILL (5, 1.0, BLC(3)) CALL RFILL (5, 1.0, TRC(3)) CALL FILL (5, 1, IDEPTH) BLC(1) = XMIN * 1000.0 TRC(1) = XMAX * 1000.0 BLC(2) = YMIN * 1000.0 TRC(2) = YMAX * 1000.0 XYRATO = (TRC(2) - BLC(2)) / (TRC(1) - BLC(1)) C Initialize plot file line drw. CALL GINITL (BLC, TRC, XYRATO, CH, IDEPTH, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GLTYPE (1, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Draw borders. CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (TRC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (TRC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C tick marks: Y axis IF (LTYPE.GT.1) THEN IF (PFLAG) GO TO 75 TICLEN = (TRC(1) - BLC(1)) / 70.0 YTICEL = BLC(1) + TICLEN YTICER = TRC(1) - TICLEN C Find interval value. DIST = FAC * (LOCRAN(2) - LOCRAN(1)) XINT = 9.0 DO 20 I = 1,21 XNOINT = AINT (DIST/XINTER(I)) IF (XNOINT.LE.XINT) GO TO 30 20 CONTINUE I = 0 C Interval and no of inter found. 30 IF (I.GT.0) THEN XINT = XINTER(I) * 1000.0 INOINT = XNOINT + 2.5 XVAL = AINT (FAC*BLC(2)/XINT) * XINT IF (XVAL.GT.FAC*BLC(2)) XVAL = XVAL - XINT IXO = I DCXM = -0.5 C Loop for all tics. DO 40 I = 1,INOINT XVAL = XVAL + XINT YPOS = XVAL / FAC IF (YPOS.LE.TRC(2)) THEN C TOP tic. CALL GPOS (TRC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (YTICER, YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Left hand tic. CALL GPOS (YTICEL, YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (BLC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C Write value. IF (MOD(ABS(LABEL),100).GT.2) THEN WRITE (MSGBUF,1030) XVAL/1000.0 CALL CHTRIM (MSGBUF, 14, MSGBUF, INCHAR) IF (IXO.GT.3) INCHAR = INCHAR - 2 DCX = - INCHAR - 1.0 DCY = -0.5 DCXM = MIN (DCXM, DCX) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, * IRET) IF (IRET.NE.0) GO TO 999 END IF END IF 40 CONTINUE END IF C X axis YTICEL = BLC(2) + TICLEN YTICER = TRC(2) - TICLEN C Determine label range DIST = FAC * (RANGE(2) - RANGE(1)) ODIST = FAC * RANGE(1) C Get interval DO 50 ITRY = 1,21 XNOINT = AINT (DIST/XINTER(ITRY)) IF (XNOINT.LE.9.0) GO TO 60 50 CONTINUE ITRY = 0 C Bottom (value) tics 60 IF (ITRY.GT.0) THEN XINT = XINTER(ITRY) * 1000.0 DCY = -1.5 INOINT = XNOINT + 2.5 XVAL = AINT (FAC*BLC(1)/XINT) * XINT IF (XVAL.GE.FAC*BLC(1)) XVAL = XVAL - XINT DO 70 I = 1,INOINT XVAL = XVAL + XINT X = XVAL / FAC IF (X.LE.TRC(1)) THEN CALL GPOS (X, YTICER, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (X, TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GPOS (X, YTICEL, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (X, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (MOD(ABS(LABEL),100).GT.2) THEN WRITE (MSGBUF,1030) XVAL/1000.0 CALL CHTRIM (MSGBUF, 14, MSGBUF, INCHAR) IF (ITRY.GT.3) INCHAR = INCHAR - 2 DCX = 0.5 - INCHAR CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, * IRET) IF (IRET.NE.0) GO TO 999 END IF END IF 70 CONTINUE END IF C Y axis label 75 DCX = DCXM - 2.0 YPOS = (TRC(2) + BLC(2)) / 2.0 CALL GPOS (BLC(1), YPOS, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (K.EQ.1) THEN MSGBUF = 'Imaginary part in ' // PREFIX // ' Jy' ELSE MSGBUF = 'Imaginary part in ' // PREFIX // ' gain' END IF CALL REFRMT (MSGBUF, '_', INCHAR) DCY = INCHAR / 2.0 - 1.0 CALL GCHAR (INCHAR, 1, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C X axis label DCY = -1.5 IF (MOD(ABS(LABEL),100).GT.2) DCY = -2.833 X = (TRC(1) + BLC(1)) / 2.0 CALL GPOS (X, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (K.EQ.1) THEN MSGBUF = 'Real part in ' // PREFIX // ' Jy' ELSE MSGBUF = 'Real part in ' // PREFIX // ' gain - 1' END IF CALL REFRMT (MSGBUF, '_', INCHAR) DCX = 0.5 - INCHAR / 2.0 CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (MOD(ABS(LABEL),100).LT.7) THEN NXFR = 0 NXST = 0 NAX = CATBLK(KIDIM) INC = 2 DO 80 I = 1,NAX ICPNT = KHCTP+(I-1)*INC CALL H2CHR (8, 1, CATH(ICPNT), CTEMP) IF (CTEMP(1:4).EQ.'FREQ') NXFR = I IF (CTEMP(1:4).EQ.'STOK') NXST = I 80 CONTINUE C Source name, stokes, freq. CALL GPOS (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 DCX = 0.0 DCY = 0.5 IANGL = 0 CALL H2CHR (8, 1, CATH(KHOBJ), CTEMP) FREQ = 0.0 JSTOK = 12 IF (NXFR.GT.2) FREQ = CATD(KDCRV+NXFR-1) + * CATR(KRCIC+NXFR-1) * (BLC(NXFR) - CATR(KRCRP+NXFR-1)) FREQ = FREQ / 1.E6 IF (NXST.GT.2) JSTOK = CATD(KDCRV+NXST-1) + 2.5 + * CATR(KRCIC+NXST-1)* (BLC(NXST) - CATR(KRCRP+NXST-1)) IF (NXFR.GT.2) WRITE (MSGBUF,1080) CTEMP, CSTOK(JSTOK), * FREQ IF (NXFR.LE.2) WRITE (MSGBUF,1080) CTEMP, CSTOK(JSTOK) CALL REFRMT (MSGBUF, '_', INCHAR) C image name INCHAR = INCHAR + 1 IF (INCHAR.GT.1) THEN MSGBUF(INCHAR:INCHAR+2) = ' _' INCHAR = INCHAR + 3 END IF CALL H2CHR (12, KHIMNO, CATH(KHIMN), NAMSTR(1:12)) CALL H2CHR (6, KHIMCO, CATH(KHIMC), NAMSTR(13:18)) CALL NAMEST (NAMSTR, CATBLK(KIIMS), MSGBUF(INCHAR:), ITMP) CALL REFRMT (MSGBUF, '_', INCHAR) CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 C time/date, version IF (LABEL.GT.0) THEN CALL GPOS (BLC(1), TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL ZDATE (ID) CALL ZTIME (IT) CALL TIMDAT (IT, ID, TIME, DATE) WRITE (MSGBUF,1085) IVER, DATE, TIME CALL REFRMT (MSGBUF, '_', INCHAR) DCY = DCY + 1.333 CALL GCHAR (INCHAR, 0, DCX, DCY, MSGBUF, BUFF1, IRET) END IF END IF END IF CALL GLTYPE (4, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 IF (K.EQ.1) THEN IF ((XMIN.LT.0.) .AND. (XMAX.GT.0.0)) THEN CALL GPOS (0.0, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (0.0, TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF ELSE IF ((XMIN.LT.0.0) .AND. (XMAX.GT.0.0)) THEN CALL GPOS (0.0, BLC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (0.0, TRC(2), BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF END IF IF ((YMIN.LT.0.0) .AND. (YMAX.GT.0.0)) THEN CALL GPOS (BLC(1), 0.0, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GVEC (TRC(1), 0.0, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 END IF CALL GLTYPE (2, BUFF1, IRET) IF (IRET.NE.0) GO TO 999 CALL GPOS (BLC(1), BLC(2), BUFF1, IRET) CALL GVEC (BLC(1), BLC(2), BUFF1, IRET) C 999 RETURN C----------------------------------------------------------------------- 1030 FORMAT (F14.1) 1080 FORMAT (A,' _',A4,'_ ',F10.3,' MHz') 1085 FORMAT ('Plot file version',I4,'__created ',A,A) END SUBROUTINE IMGDRW (NP, MAP, NSM, S1M, S2M, XMIN, XMAX, YMIN, YMAX, * IGBLK, IERR) C----------------------------------------------------------------------- C CONDRW will write commands to a plot file for the execution of C a contour plot. C Inputs: C NP I Number points on a side C MAP D(*) Image C NSM I Size of smoothing kernel C IGBLK I(256) I/O block for graph file. C Common: C CATBLK I(256) map header. C CNTRBU R(8192) buffers C Output: C IERR I error code. 0 = ok. C----------------------------------------------------------------------- INTEGER NP, NSM, IGBLK(256), IERR DOUBLE PRECISION MAP(NP,*) REAL XMIN, XMAX, YMIN, YMAX, S1M(NP,*), S2M(NP,*) C REAL VAL(3), XPOS(3), YPOS(3), TEMP, VC, VL, VM, VS, XA, XB, * XL, XLAST, XM, XS, YA, YB, YL, YLAST, YM, YS, DELTAX, DELTAY, * TLEV, XLEV, ALEVS(30), XP, YP, XSCALE, YSCALE INTEGER IPERM(3,6), IBLCX, IBLCY, IBLCY1, ICOL, II, INDEX, * INLEVS, INPIXS, IPLUS, IPOS, IROW, ISKIP, ITRCX, ITRCXM, * ITRCY, ITRI, I, MININT, LOCINT, IROUND, ISLEV, JJ, LASTC DOUBLE PRECISION DMAX INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DDCH.INC' INCLUDE 'INCS:DCAT.INC' INCLUDE 'INCS:DCNT.INC' INCLUDE 'INCS:DGPH.INC' INCLUDE 'INCS:DMSG.INC' DATA IPERM /1,3,2, 3,2,1, 3,1,2, 2,1,3, 1,2,3, 2,3,1/ C----------------------------------------------------------------------- C smooth input CALL IMGSMT (NSM, NP, MAP, S1M, S2M) XSCALE = (XMAX - XMIN) / (NP -1.0) YSCALE = (YMAX - YMIN) / (NP -1.0) ISKIP = 2 XLAST = -1000. YLAST = -1000. TEMP = NP - 1 TEMP = 10 - 3 * LOG10 (TEMP) LOCINT = IROUND (TEMP) IF (LOCINT.LT.2) LOCINT = 2 IBLCY = 1 ITRCY = NP ITRCX = NP IBLCX = 1 INPIXS = ITRCX - IBLCX + 1 LASTC = 0 C magic parms for dashed lines XLEV = 256.0 / INPIXS TLEV = 256.0 / (ITRCY - IBLCY + 1.0) ISLEV = SQRT (1.0 / (XLEV * TLEV)) + 0.1 IF (ISLEV.LT.1) ISLEV = 1 IF (XLEV.LT.1.0) XLEV = (SQRT (XLEV) + 3.0*XLEV) / 4.0 C Determine number of levels. C and convert to absolute levels. DMAX = 0.0D0 DO 20 IROW = 1,NP DO 10 I = 1,NP DMAX = MAX (DMAX, MAP(I,IROW)) 10 CONTINUE 20 CONTINUE INLEVS = 2.0D0 * LOG10 (DMAX) DO 45 I = 1,30 ALEVS(I) = I/2.0 45 CONTINUE WRITE (MSGTXT,1045) DMAX, ALEVS(INLEVS) CALL MSGWRT (3) C Read and save first row. DO 55 I = IBLCX,ITRCX IF (MAP(I,IBLCY).GT.0.0D0) THEN RLROW(I-IBLCX+1) = LOG10 (MAP(I,IBLCY)) ELSE RLROW(I-IBLCX+1) = -0.3 END IF 55 CONTINUE C Loop over all rows. IBLCY1 = IBLCY + 1 DO 300 IROW = IBLCY1,ITRCY C Read proper row. DO 110 I = IBLCX,ITRCX IF (MAP(I,IROW).GT.0.0D0) THEN BUFF(I) = LOG10 (MAP(I,IROW)) ELSE BUFF(I) = -0.3 END IF 110 CONTINUE IPOS = IBLCX C Loop over all pixels in row. IPLUS = 0 ITRCXM = ITRCX - 1 DO 250 ICOL = IBLCX,ITRCXM IPLUS = IPLUS + 1 C Init values VAL(1) = BUFF(IPOS+IPLUS-1) VAL(2) = BUFF(IPOS+IPLUS) VAL(3) = RLROW(IPLUS) C Init positions. XPOS(1) = ICOL XPOS(2) = ICOL + 1 XPOS(3) = ICOL YPOS(1) = IROW YPOS(2) = IROW YPOS(3) = IROW - 1 C Loop for both triangles. DO 200 ITRI = 1,2 C Changes for 2nd triangle. IF (ITRI.EQ.2) THEN VAL(1) = RLROW(IPLUS+1) XPOS(1) = ICOL + 1 YPOS(1) = IROW - 1 END IF C Order points in triangle. DO 130 II = 1,3 IF (VAL(II).EQ.FBLANK) GO TO 200 130 CONTINUE INDEX = 0 IF (VAL(1).GT.VAL(2)) INDEX = 1 IF (VAL(3).GE.VAL(1)) INDEX = INDEX + 2 IF (VAL(2).GE.VAL(3)) INDEX = INDEX + 4 C find large, med, small C values and X,Y positions. II = IPERM(1,INDEX) VL = VAL(II) XL = XPOS(II) YL = YPOS(II) C II = IPERM(2,INDEX) VM = VAL(II) XM = XPOS(II) YM = YPOS(II) C II = IPERM(3,INDEX) VS = VAL(II) XS = XPOS(II) YS = YPOS(II) C Loop for all levels. DO 190 II = 1,INLEVS VC = ALEVS(II) C Cut down negatives IF (VC.GE.0.0) GO TO 140 IF ((XLEV.LT.2.85) .AND. (ITRI.EQ.2)) GO TO 190 IF (XLEV.GE.1.0) GO TO 140 JJ = IROW + ICOL + II IF (MOD(JJ, ISLEV).NE.0) GO TO 190 140 IF ((VC.GT.VL) .OR. ((VL-VS).LE.0.0)) GO TO 200 C If level not right, next lev. IF (VC.LE.VS) GO TO 190 C Interpolate btwn max two corns. TEMP = (VC-VS) / (VL-VS) XA = TEMP * (XL-XS) + XS YA = TEMP * (YL-YS) + YS C See which corners 2nd pt. btwn. IF (VC.GT.VM) GO TO 150 IF (VM.EQ.VS) GO TO 150 C Level btwn med & small corners. TEMP = (VC-VS) / (VM-VS) XB = TEMP * (XM-XS) + XS YB = TEMP * (YM-YS) + YS GO TO 160 C Level btwn large & med corners. 150 TEMP = (VC-VM) / (VL-VM) XB = TEMP * (XL-XM) + XM YB = TEMP * (YL-YM) + YM C Issue position & write commands C We can avoid position command C if we switch A and B. 160 IF ((XLAST.EQ.XB) .AND. (YLAST.EQ.YB)) THEN TEMP = XA XA = XB XB = TEMP TEMP = YA YA = YB YB = TEMP C See if we need to position. ELSE IF ((XLAST.NE.XA) .OR. (YLAST.NE.YA)) THEN XP = (XMIN + (XA - 1.) * XSCALE) * 1000.0 YP = (YMIN + (YA - 1.) * YSCALE) * 1000.0 CALL GPOS (XP, YP, IGBLK, IERR) IF (IERR.NE.0) GO TO 999 END IF C Draw vector. IF (VC.GE.0.0) THEN XP = (XMIN + (XB - 1.) * XSCALE) * 1000.0 YP = (YMIN + (YB - 1.) * YSCALE) * 1000.0 CALL GVEC (XP, YP, IGBLK, IERR) IF (IERR.NE.0) GO TO 999 XLAST = XB YLAST = YB GO TO 190 C Negative contours broken ELSE TEMP = LOCINT * SQRT (((XB-XA)**2 + (YB-YA)**2) * / 2.0) MININT = IROUND (TEMP) IF (MININT.LT.2) MININT = 2 DELTAX = (XB - XA) / MININT DELTAY = (YB - YA) / MININT DO 185 I = 1,MININT,2 XB = XA + DELTAX YB = YA + DELTAY XP = (XMIN + (XB - 1.) * XSCALE) * 1000.0 YP = (YMIN + (YB - 1.) * YSCALE) * 1000.0 CALL GVEC (XP, YP, IGBLK, IERR) IF (IERR.NE.0) GO TO 999 IF (I.LT.MININT-1) THEN XA = XB + DELTAX YA = YB + DELTAY XP = (XMIN + (XA - 1.) * XSCALE) * 1000.0 YP = (YMIN + (YA - 1.) * YSCALE) * 1000.0 CALL GPOS (XP, YP, IGBLK, IERR) IF (IERR.NE.0) GO TO 999 END IF 185 CONTINUE XLAST = XB YLAST = YB END IF 190 CONTINUE 200 CONTINUE 250 CONTINUE CALL RCOPY (INPIXS, BUFF(IPOS), RLROW) 300 CONTINUE C 999 RETURN C----------------------------------------------------------------------- 1045 FORMAT ('Gridded image max=',1PE12.4,' counts; peak contour', * 0PF4.1,' in log10') END SUBROUTINE IMGSMT (NSM, NP, MAP, SUM, WGT) C----------------------------------------------------------------------- C Inputs: C NP I size of image C In/out: C MAP D(*) image: in unsmoothed, out smoothed C Output: C SUM R(*) Sum C WGT R(*) Weight sum C----------------------------------------------------------------------- INTEGER NP, NSM DOUBLE PRECISION MAP(NP,*) REAL SUM(NP,*), WGT(NP,*) C INTEGER NK, NC INTEGER I, J, IX, IY, N REAL KERNEL(33,33), W, R C----------------------------------------------------------------------- IF (NSM.LE.2) GO TO 999 N = NP * NP CALL RFILL (N, 0.0, SUM) CALL RFILL (N, 0.0, WGT) NK = MIN (33, NSM) NK = (NK / 2) * 2 + 1 NC = NK / 2 + 1 W = NK / 2 C kernel DO 20 I = 1,NK DO 10 J = 1,NK R = (I-NC)*(I-NC) + (J-NC)*(J-NC) KERNEL(J,I) = EXP (-R/W) 10 CONTINUE 20 CONTINUE C sum DO 90 IY = 1,NP DO 80 IX = 1,NP DO 40 I = MAX(1,IY-NC+1),MIN(NP,IY+NC-1) DO 30 J = MAX(1,IX-NC+1),MIN(NP,IX+NC-1) W = KERNEL(J-IX+NC,I-IY+NC) SUM(J,I) = SUM(J,I) + W * MAP(IX,IY) WGT(J,I) = WGT(J,I) + W 30 CONTINUE 40 CONTINUE 80 CONTINUE 90 CONTINUE C average DO 120 IY = 1,NP DO 110 IX = 1,NP IF (WGT(IX,IY).GT.0.0) THEN MAP(IX,IY) = SUM(IX,IY) / WGT(IX,IY) ELSE MAP(IX,IY) = 0.0D0 END IF 110 CONTINUE 120 CONTINUE C 999 RETURN END SUBROUTINE EVAUVH C----------------------------------------------------------------------- C EVAUVH writes HI files and copies extension files for the C subtracted and the divided data sets when they are kept C----------------------------------------------------------------------- C INTEGER I, J, IERR, DISK, CNO, CATB(256), LUN1, LUN2, NONOT CHARACTER NAMO*12, CLAO*6, AOPCOD(2)*4, HILINE*72, AMETH*4, * AMODL*4, NOTTYP*2 HOLLERITH CATBH(256) INCLUDE 'EVAUV.INC' INCLUDE 'INCS:DMSG.INC' INCLUDE 'INCS:DFIL.INC' INCLUDE 'INCS:DGDS.INC' INCLUDE 'INCS:DHDR.INC' INCLUDE 'INCS:DCAT.INC' EQUIVALENCE (CATB, CATBH) DATA LUN1, LUN2 /27,28/ DATA AOPCOD /'SUB','DIV'/ DATA NONOT, NOTTYP /0,' '/ C----------------------------------------------------------------------- CALL HIINIT (3) DISK = FVOL(-SC1FIL) CNO = FCNO(-SC1FIL) DO 100 I = 1,2 C get image header CALL CATIO ('READ', DISK, CNO, CATBLK, 'REST', BUFF1, IERR) IF ((IERR.GT.0) .AND. (IERR.LT.5)) THEN WRITE (MSGTXT,1000) IERR, 'READING HEADER', I GO TO 70 END IF CALL HISCOP (LUN1, LUN2, DISKIN, DISK, OLDCNO, CNO, CATBLK, * BUFF1, BUFF2, IERR) IF (IERR.GT.2) THEN WRITE (MSGTXT,1000) IERR, 'COPYING HI FILE' GO TO 70 END IF C calibration history CALL CALHIS (LUN2, BUFF2, IERR) IF (IERR.NE.0) GO TO 100 C New history CALL HENCO1 (TSKNAM, NAMEIN, CLAIN, SEQIN, DISKIN, LUN2, * BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C If point model, no model file IF (.NOT.DOPTMD) THEN C Model file(s) CALL HENCO2 (TSKNAM, NAME2, CLAS2, SEQ2, DISK2, LUN2, BUFF2, * IERR) IF (IERR.NE.0) GO TO 75 WRITE (HILINE,1020) TSKNAM, VER IF (MODEL.EQ.1) CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C Number of input images WRITE (HILINE,1025) TSKNAM, MFIELD CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C Add no. clean comps. DO 25 J = 1,MFIELD NCOMP(J) = NSUBG(J) - 1 WRITE (HILINE,1021) TSKNAM, J, BCOMP(J), J, NCOMP(J) IF (MODEL.EQ.1) CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 25 CONTINUE END IF CALL H2CHR (12, KHIMNO, CATH(KHIMN), NAMO) CALL H2CHR (6, KHIMCO, CATH(KHIMC), CLAO) CALL HENCOO (TSKNAM, NAMO, CLAO, CATB(KIIMS), DISK, LUN2, * BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C Modeling method AMETH = 'DFT ' IF (METHOD.EQ.1) AMETH = 'GRID' WRITE (HILINE,1027) TSKNAM, AMETH CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C OPCODE WRITE (HILINE,1029) TSKNAM, AOPCOD(I) CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C Model type C CC or image IF (.NOT.DOPTMD) THEN AMODL = 'COMP' IF (MODEL.EQ.2) AMODL = 'IMAG' IF (MODEL.EQ.3) AMODL = 'SUBI' WRITE (HILINE,1028) TSKNAM, AMODL CALL HIADD (LUN2, HILINE, BUFF2, IERR) IF (IERR.NE.0) GO TO 75 C Point model ELSE WRITE (HILINE,1030) TSKNAM, SMODEL(1), SMODEL(2), SMODEL(3) CALL HIADD (LUN2, HILINE, BUFF3, IERR) IF (IERR.NE.0) GO TO 75 C Other parameters WRITE (HILINE,1031) TSKNAM, SMODEL(4), SMODEL(5), SMODEL(6), * SMODEL(7) IF (SMODEL(4).GT.0.01) CALL HIADD (LUN2, HILINE, BUFF3, * IERR) IF (IERR.NE.0) GO TO 75 END IF GO TO 75 70 CALL MSGWRT (7) C Close history file. 75 CALL HICLOS (LUN2, .TRUE., BUFF2, IERR) C Copy tables C CALL ALLTAB (NONOT, NOTTYP, LUN1, LUN2, DISKIN, DISK, OLDCNO, C * CNO, CATBLK, BUFF1, BUFF2, IERR) CALL COPTAB (DISKIN, OLDCNO, DISK, CNO, IERR) IF (IERR.GT.2) THEN WRITE (MSGTXT,1000) IERR, 'COPYING TABLES', I CALL MSGWRT (7) END IF C Update CATBLK. CALL CATIO ('UPDT', DISK, CNO, CATBLK, 'REST', BUFF1, IERR) DISK = FVOL(-SC2FIL) CNO = FCNO(-SC2FIL) 100 CONTINUE C 999 RETURN C----------------------------------------------------------------------- 1000 FORMAT ('EVAUVH: ERROR',I4,' ON ',A,' OUTPUT FILE',I2) 1020 FORMAT (A6,'VER = ',I6,' / CC file ver. no.') 1021 FORMAT (A6,'BCOMP(',I2,')=',I6,', NCOMP(',I3,') =',I6, * ' / First-last comp. no.') 1025 FORMAT (A6,'NMAPS = ',I6,' / Number of model images') 1027 FORMAT (A6,'CMETHOD = ''',A4,''' / Model method') 1028 FORMAT (A6,'CMODEL = ''',A4,''' / Model type') 1029 FORMAT (A6,'OPCODE = ''',A4,''' / Operation type') 1030 FORMAT (A6,'SMODEL = ',F12.5,2F10.5,' / Model flux,RA,Dec') 1031 FORMAT (A6,' ',F12.5,2F10.5,F8.5,' / Other parms') END