SUBROUTINE VSAA * * Module number: * * Module name: * * Keyphrase: * ---------- * Process particle monitor data obtained from HSP * * Description: * ------------ * Read particle monitor data (dark counts) from HSP and put the accumulated * data in three output files: (1) number of observations (2) accumulated * count rate (3) accumulated count rate square * input files should be single-group files. * * FORTRAN name: VSAA.FOR * * Keywords of accessed files and tables: * -------------------------------------- * Name I/O Description / Comments * * input parameters * * 'infiles' I input file template name * 'inmasks' I input mask template name * 'nfiles' I number of input files * 'intable' I input state vector table name * 'outfile1', 'outfile2', 'outfile3' * I output file(s) containing accumulated * result * 'npoints' I number of points used in Lagrange's * polynomial interpolation * * input parameters from file header * * 'SAMPTIME' I sample time * 'FPKTTIME' I start time of the observation * 'LPKTTIME' I stop time of the observation * 'CD1_1' I interval between each data point * 'DQMSKFV' I bad pixel mask value * 'DATA_TYP' I data type * 'TRUE_CNT' I digital data calibration switch * * Subroutines Called: * ------------------- * CDBS: * VGTXYZ, VGTSAA, VXYZIN, VMJD * SDAS: * UCLGSI, UCLGSR, UCLGST, UIMOPN, UIMGID, UIMCLO, UTRGTR, * UIGL2R, UIPL2R, UTTCLO, UMSPUT, UHDGSR, UHDGSD * Others: * None * * History: * -------- * Version Date Author Description * 1 03-22-89 J.-C. HSU design and coding *------------------------------------------------------------------------------- * *== local: * --total number of valid input data points INTEGER NPTS, * --error status : STATUS, * --image/mask descriptor : IMGDSC, MSKDSC, * --number of dimensions : NAXIS, NAX, * --maximum size of the state vector table : TSIZE, * --size of the output files : NX, NY, * --array sizes : DIMEN(7), DIM(7), * --number of input files : NFL, NFILES, * --input file/mask ID : INFID, INMID, * --data type : DTYPE, TYPE, * --dynamic memory pointer of input file * --and mask : FILEPT, MASKPT, * --image descripters and dynamic memory * --pointer of output files : OUTID(3), OUTPT(3), * --number of points used in Lagrange's * --polynomial interpolation : NLAG, * --loop indices : I, J, K, IX, IY, INDX, * --status : STAT(20), STATOK * --data and mask buffer PARAMETER (TSIZE = 2000) PARAMETER (NX = 360) PARAMETER (NY = 61) * --time between data points REAL INTRVL, * --mask fill values : FILVAL, * --integration time : SMPTM, * --longitude and latitude : LONG, LAT, * --count rate : CTRATE * --epoch in the state vector table DOUBLE PRECISION MJD(TSIZE), * --Greenwich mean sidereal time, right * --ascension, declination, and geocentric * --distance : GMST, RA, DECL, R, * --XYZ positions of the state vectors : XYZ(3, TSIZE), PXYZ0(3), PXYZ(3), : EPOCH, STRTJD, STOPJD, MJ2000 * --calibration switch CHARACTER*8 CALIB * --data type CHARACTER*10 DAT * --time tag CHARACTER*24 STRTTM, STOPTM * --input file, mask, table, output file, * --input file template names CHARACTER*128 IFILE, ITABLE, IMASK, FNAMES(3), : INFILE, INMASK * --error message context CHARACTER*130 CONTXT, MESS PARAMETER (MJ2000 = 51544.5D0) *==========================begin iraf77.inc (without INTEGER*2)================= * Include file for the iraf77 FORTRAN interface to the IRAF VOS * Get IRAF common into main program * LOGICAL MEMB(1) INTEGER MEMI(1) INTEGER MEML(1) REAL MEMR(1) DOUBLE PRECISION MEMD(1) COMPLEX MEMX(1) EQUIVALENCE (MEMB, MEMI, MEML, MEMR, MEMD, MEMX) COMMON /MEM/ MEMD * * File I/O access modes * INTEGER RDONLY PARAMETER (RDONLY = 1) INTEGER RDWRIT PARAMETER (RDWRIT = 2) INTEGER WRONLY PARAMETER (WRONLY = 3) INTEGER APPEND PARAMETER (APPEND = 4) INTEGER NEWFIL PARAMETER (NEWFIL = 5) INTEGER TMPFIL PARAMETER (TMPFIL = 6) INTEGER NEWCPY PARAMETER (NEWCPY = 7) INTEGER NEWIMG PARAMETER (NEWIMG = 5) INTEGER EOF PARAMETER (EOF = -2) * * codes for data types * INTEGER TYBOOL PARAMETER (TYBOOL = 1) INTEGER TYCHAR PARAMETER (TYCHAR = 2) INTEGER TYSHOR PARAMETER (TYSHOR = 3) INTEGER TYINT PARAMETER (TYINT = 4) INTEGER TYLONG PARAMETER (TYLONG = 5) INTEGER TYREAL PARAMETER (TYREAL = 6) INTEGER TYDOUB PARAMETER (TYDOUB = 7) INTEGER TYCPLX PARAMETER (TYCPLX = 8) INTEGER TYUSHT PARAMETER (TYUSHT = 11) INTEGER TYUBYT PARAMETER (TYUBYT = 12) * * TYTEXT is a special code for the iraf77 interface; it is not in the VOS * INTEGER TYTEXT PARAMETER (TYTEXT = 13) *========================end iraf77.inc========================================= *=========================begin hsp.inc========================================= * --status return code INTEGER OK, ERRNUM(20) * --message destination and priority INTEGER DEST, PRIO DATA OK /0/ DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710, : 711, 712, 713, 714, 715, 716, 717, 718, 719, 720/ DATA DEST, PRIO /1, 0/ *=========================end hsp.inc=========================================== *------------------------------------------------------------------------------ * * read parameter file * CALL UCLGST ('infiles', INFILE, STAT(1)) CALL UCLGST ('inmasks', INMASK, STAT(2)) CALL UCLGSI ('nfiles', NFILES, STAT(3)) CALL UCLGST ('intable', ITABLE, STAT(4)) CALL UCLGST ('outfile1', FNAMES(1), STAT(5)) CALL UCLGST ('outfile2', FNAMES(2), STAT(6)) CALL UCLGST ('outfile3', FNAMES(3), STAT(7)) CALL UCLGSI ('npoints', NLAG, STAT(8)) * DO 10 I = 1, 8 IF (STAT(I) .NE. OK) THEN STATUS = ERRNUM(1) CONTXT = 'can not get input parameter(s)' GO TO 999 END IF 10 CONTINUE * * read the input state vector table * CALL VGTXYZ (ITABLE, TSIZE, MJD, XYZ, NPTS, STATUS) * IF (STATUS .NE. OK) THEN CONTXT = 'error reading input state vector table' GO TO 999 END IF * * read/create output (accumulation) files * CALL VGTSAA (FNAMES, 3, NX, NY, OUTID, OUTPT, STATUS) * IF (STATUS .NE. OK) THEN CONTXT = 'error reading output file(s)' GO TO 999 END IF * * initialize template processing for both input data files and masks * CALL UIMOTP (INFILE, IMGDSC, STAT(1)) CALL UIMOTP (INMASK, MSKDSC, STAT(2)) IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN STATUS = ERRNUM(3) CONTXT = 'cannot initialize template processing' GO TO 999 END IF * * process each input file * NFL = 0 DO 100 I = 1, NFILES * * expand image name template * CALL UIMXTP (IMGDSC, IFILE, STAT(1)) CALL UIMXTP (MSKDSC, IMASK, STAT(2)) * * if end of file list, exit from the loop * IF (STAT(1) .EQ. EOF .OR. STAT(2) .EQ. EOF) GO TO 200 * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN STATUS = ERRNUM(1) CONTXT = 'cannot expand file name template' GO TO 999 END IF * * open input file * CALL UIMOPN (IFILE, RDONLY, INFID, STAT(1)) CALL UIMGID (INFID, TYPE, NAX, DIM, STAT(2)) * * if cannot open file, skip, and process the next file * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN CONTXT = 'VSAA: cannot open input file or get image ' : // 'description' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 100 END IF * * open input mask file * CALL UIMOPN (IMASK, RDONLY, INMID, STAT(1)) CALL UIMGID (INMID, DTYPE, NAXIS, DIMEN, STAT(2)) * * if cannot open mask, skip, and process the next file * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN CONTXT = 'VSAA: cannot open input mask or get image ' : // 'description' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 100 END IF * * check mask image description against those of the data file * IF (DTYPE .NE. TYPE .OR. NAXIS .NE. NAX .OR. : DIMEN(1) .NE. DIM(1) ) THEN CONTXT = 'VSAA: input mask has ' : // 'incompatible data type or dimension or size' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF * * read keywords from input data file header * CALL UHDGSR (INFID, 'SAMPTIME', SMPTM, STAT(1)) CALL UHDGST (INFID, 'FPKTTIME', STRTTM, STAT(2)) CALL UHDGST (INFID, 'LPKTTIME', STOPTM, STAT(3)) CALL UHDGSR (INFID, 'CD1_1', INTRVL, STAT(4)) CALL UHDGSR (INFID, 'DQMSKFV', FILVAL, STAT(5)) CALL UHDGST (INFID, 'DATA_TYP', DAT, STAT(6)) CALL UHDGST (INFID, 'TRUE_CNT', CALIB, STAT(7)) * DO 20 J = 1, 7 IF (STAT(J) .NE. OK) THEN CONTXT = 'VSAA: error reading header keyword(s)' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF 20 CONTINUE * * translate date/time string to MJD * CALL VMJD (STRTTM, STRTJD, STAT(1)) CALL VMJD (STOPTM, STOPJD, STAT(2)) * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN STATUS = ERRNUM(1) CONTXT = 'VSAA: error decoding time string' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF * * check if the file epoch is within range of the state vector table * IF (STRTJD .LT. MJD(1) .OR. STRTJD .GT. MJD(NPTS) .OR. : STOPJD .LT. MJD(1) .OR. STOPJD .GT. MJD(NPTS)) THEN STATUS = ERRNUM(1) CONTXT = 'VSAA: file epoch outside state vector table ' : // 'range' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF * * allocate dynamic memory for the input data * CALL UDMGET (DIM(1), TYREAL, FILEPT, STAT(1)) CALL UDMGET (DIM(1), TYREAL, MASKPT, STAT(2)) * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN CONTXT = 'VSAA: error allocating dynamic memory' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF * * read data from input file and its mask * CALL UIGL1R (INFID, MEMR(FILEPT), STAT(1)) CALL UIGL1R (INMID, MEMR(MASKPT), STAT(2)) * IF (STAT(1) .NE. OK .OR. STAT(2) .NE. OK) THEN CONTXT = 'VSAA: error reading data from input file ' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 40 END IF * * process each point * DO 30 K = 1, DIM(1) * * skip masked points * IF (MEMR(MASKPT+K-1) .NE. FILVAL) THEN * * figure out the epoch of the data point * EPOCH = STRTJD + ((K-1) * INTRVL + 0.5 * SMPTM) / : 86400.D0 * * use Lagrangian interpolation to find the state vector corresponding to the * epoch of the data point * CALL VXYZIN (MJD, XYZ, NPTS, NLAG, EPOCH, PXYZ0, : STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'VSAA: error in Lagrangian polynomial ' : // 'interpolation' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) GO TO 50 END IF * * convert the state vector XYZ to longitude and latitude * * precess the state vector from J2000 to the current epoch * CALL CDPRC (PXYZ0, MJ2000, EPOCH, PXYZ) * * convert XYZ to right ascension and declination (=latitude) in decimal degrees * CALL CDPCOR (PXYZ, R, RA, DECL) LAT = REAL (DECL) * * find the mean Greenwich Mean Sidereal Time (in degrees) * CALL CDLMST(EPOCH, 0., GMST) GMST = GMST * 15.D0 * * calculate logitude and put within the range of [0, 360) * LONG = REAL (RA - GMST) LONG = MOD (LONG, 360.) IF (LONG .LT. 0.) LONG = LONG + 360. * * determine the pixel ID from longitude and latitude * IX = INT(LONG) + 1 IY = INT(LAT) + 31 INDX = 360 * (IY - 1) + IX * * data type is digital, divide by integration time * CTRATE = MEMR(FILEPT+K-1) / SMPTM MEMR(OUTPT(1)+INDX-1) = MEMR(OUTPT(1)+INDX-1) +1 MEMR(OUTPT(2)+INDX-1) = MEMR(OUTPT(2)+INDX-1) + : CTRATE MEMR(OUTPT(3)+INDX-1) = MEMR(OUTPT(3)+INDX-1) + : CTRATE ** 2 END IF 30 CONTINUE * * release dynamic memory * 40 CALL UDMFRE (FILEPT, TYREAL, STAT(1)) CALL UDMFRE (MASKPT, TYREAL, STAT(2)) * * close input file/mask * 50 CALL UIMCLO (INFID, STAT(3)) CALL UIMCLO (INMID, STAT(4)) DO 60 J = 1, 4 IF (STAT(J) .NE. OK) THEN CONTXT = 'cannot free dynamic memory or close ' : // 'input file ' CALL UMSPUT (CONTXT, DEST, PRIO, STATOK) END IF 60 CONTINUE * NFL = NFL + 1 100 CONTINUE * * close image list * 200 CALL UIMCTP (IMGDSC, STATOK) CALL UIMCTP (MSKDSC, STATOK) * IF (NFL .LE. 0) THEN STATUS = ERRNUM(2) CONTXT = 'no input data been successfully read' GO TO 999 END IF * * update output files * * * write the result to output file/mask * DO 220 I = 1, 3 DO 210 K = 1, NY J = NX * (K-1) CALL UIPL2R (OUTID(I), K, MEMR(OUTPT(I)+J), STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot write to output file(s)' GO TO 999 END IF 210 CONTINUE 220 CONTINUE * * close output file * DO 230 J = 1, 3 CALL UIMCLO (OUTID(J), STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot close output file(s)' GO TO 999 END IF 230 CONTINUE * * deallocate output file dynamic memory * DO 240 J = 1, 3 CALL UDMFRE (OUTPT(J), TYREAL, STATUS) IF (STATUS .NE. OK) THEN CONTXT = 'cannot deallocate dynamic memory' GO TO 999 END IF 240 CONTINUE * STATUS = OK GO TO 1000 * * write error message * 999 MESS = 'VSAA: ' // CONTXT CALL UMSPUT (MESS, DEST, PRIO, STATOK) * 1000 RETURN END