idl routines to process mock ips data

Last modified: Thu Jul 22 11:48:21 2021.


List of Routines


Routine Descriptions

A_IPSEXAMPLE - LOOKING AT MOCK IPS DATA

[Next Routine] [List of Routines]
NAME:
a_ipsexample - Looking at mock ips data
   
       IPS observations use On,off position switching and the mock spectrometers
 to record data. The ips_xxxx  idl routines will input and process this data.
 The telescope/datataking setup is normally:
 - 180 second on source , 180 seconds off source.
 - The mock setup:
   327- 53 Mhz bw, 2048 channels, 2ms sampling
   lbw- 7*80Mhz bw, 2048 channels, 1ms sampling
   sbw- 4*80Mhz bw, 2048 channels, 2ms sampling
   sbh- 7*80Mhz bw, 2048 channels, 2ms sampling
   cb - 7*160Mhz bw, 4096 channels, 2ms sampling

 to process multiple on/off pairs from a single receiver for a given day:
 (this can be later expanded to multiple receivers)
 .. for this example lets do lbw
 - find a computer with lots of memory > 16GBytes
 idl71
 @phil
 @ipsinit
 yymmdd=20180824
 proj='a3285'
 freqResHz=.1   ; generate spectra with .1 hz resolution
 bwSB=10.       ; use 10Mhz sub bands

; scan all the patterns of this day, project

 n=ipsscanpatterns(projId,yymmdd,patI,verb=1)

; limit to lbw data

 ii=where(patI.sumI.h.frontend eq 'lbw',nlbw)

; 1 on/off pattern generates 7 separate files (one for each 80 Mhz band).
; use the patid to group these together

; get the unique pattern id's
  upatid=patI[uniq(patI.patid,sort(patI.patid))].patid
  npat=n_elements(upatid)
; loop over each on/off set
  firstTime=1  
  icur=0l 		; to keep track of current one we're doing
  for ipat=0,npat-1 do begin
     ii=where(patI.patid eq upatid[ipat],bandsInPat)
;    process each band
    for ib=0,bandsInPat-1 do begin
;       get the data
		 istat=ipsget(patI[ii[ib]].filelist[0],patI[ii[ib]].filelist[1],ipsI,$
                        bon=bon,boff=boff,bfitoff=bfitoff,verb=verb,bwsubband=bwSB) &$
;       compute the spectra
       ipscmpspectra,ipsI,freqResHz,ipsSpcI

;     if first time, allocate arrays to hold all the data

      if firstTime then begin
          ipsIlbwAr=replicate(ipsI,npat*bandsInPat)
          ipsspcIlbwAr=replicate(ipsspcI,npat*bandsInPat)
          firstTime=0
      endif
      ipsIlbwAr[icur]=ipsI
      ipsspcIlbwAr[icur]=ipsspcI
      icur++
    endfor  ; end bands in pattern
  endfor    ; end on/off loop

 The ipsget() routine returns an array of ipsI structs. these
Hold the total power info after the rfi removal by frquency channel
and then total power after removing narrow spikes. (see ipsget()

 The return structure contains:
; there were 3 lbw patterns of 7 bands in each pattern

 help,ipsIlbwAr
 IPSILBWAR       STRUCT    = ->  Array[21]

 help,ipsIlbwAr,/st
** Structure <282f488>, 21 tags, length=35759736, data length=35759719, refs=1:
 SRC           STRING    'B1005+077'              .. source name
 FON          STRING    '/share/pdata1/pdev/a3285.20180824.b0s1g0.01200.fits' ..onSrc
 FOFF         STRING    '/share/pdata1/pdev/a3285.20180824.b0s1g0.01300.fits' ..offSrc
 NREC         LONG               180              .. secs in on,off
 SMPREC       LONG              1000              .. 1000 samples/sec 1ms sampling
 SMPTOT       LONG            180000              .. total samples 180 secs
 BWTOT        DOUBLE           80.000000          .. bwMhz this band
 BWSUBBAND    FLOAT           10.0000             .. bwMhz for sub bands
 NSUBBAND     LONG                 8              .. # of sub bands 
 RFCFR        DOUBLE           1160.0000          .. rf centerfreqMhz for band
 RFFREQ       DOUBLE    Array[2048]               .. freq for each rf freq channel
 CFRSUBBAND   DOUBLE    Array[8]                  .. centerFreqMhz for each 10Mhz subband
 HDR          STRUCT    -> MASFHDR Array[1]       .. fits header from first on rec
 BLKLEN       INT             10                  .. blocklen (secs) for rms  rfi removal
 NBLKS        LONG                18              .. number of 10 sec blocks in 180 secs
 BLKION       STRUCT    ->  Array[18]  .. onSource totalPower info for each 10secs
 BLKIOFF      STRUCT    ->  Array[18]  .. offSource totalPower info for each 10 secs
 SPIKECLIPLEV FLOAT         4.000                 .. narrow spike clipLevel is 4 sigma 
 NCLIPON      LONG      Array[2, 8]               .. # on Narrow spikes clipped [pola/b,nsb]
 NCLIPOFF     LONG      Array[2, 8]               .. #off Narrow spikes clipped
 TPDIF        FLOAT     Array[180000, 2, 8]       .. OnSrcTp -OffSrcTp after rms rfi removal
                                                     and narrow spike removal
                                                     [TotSmpls,2pol,8subBands]

; Inside ipsIlbwAr,  the BlkIOn,BlkIOff array of structs hold the rfi removal results
; for each 10 second block.

 IDL> help,IPSILBWAR.blkion,/st
** Structure <2325778>, 6 tags, length=672836, data length=672836, refs=2:
 I0REC     INT              0       ; index first 1 sec rec used (cnt from 0)
 I1REC     INT              9       ; index last 1 sec rec used
 RMSBYCHAN FLOAT     Array[2048, 2] ; rms/mean by frequency channel [2048Channel,2pols]
 MASK      FLOAT     Array[2048, 2] ; mask for rms.1--> use, 0--> don't use 
 NGOOD     LONG      Array[2, 8]    ; # of good points each subband, pol
 TP        FLOAT     Array[10000, 2, 8]; computed total power after excluding bad freq channels
                                     10000 smpIn10secs, 2pol, 8 subbands 

;
; The ipscmpspectra, routine computes the spectra of the total power. See ipscmpspectra
;for a description of how it does it. The spectra are normalized to the source deflection.
; 
; The ipsspcI struct returned contains:

 help,ipsspcIlbwar
 IPSSPCILBWAR    STRUCT    = ->  Array[21]

IDL> help,ipsspcIlbwar,/st
** Structure <2402238>, 9 tags, length=6440088, data length=6440080, refs=1:
 SPCLEN  LONG    5000        ; length of spectra
 SPCFREQ DOUBLE  Array[5000] ; freq for each channel
 NSPC    LONG    18          ; number of 10 sec spectra
 NSB     LONG    8           ; number of 10Mhz subbands
 BWSB    FLOAT   10.0000     ; bw each subBand [MHz]
 RFCFRSB DOUBLE  Array[8]    ; rf center freq each subBand
 SPC     FLOAT   Array[5000, 18, 2, 8] ; spectra [5000 chan,18x10secs,2 pol,8 subbands]
 SPCAVG  FLOAT   Array[5000, 2, 8]     ; average the 18 10 sec spectra
 SPCMED  FLOAT   Array[5000, 2, 8]     ; take the median of the 18 10 sec spectra

; from within idl
 explain,ipsdoc   .. will list the ips routines
 explain,ipsget   .. will list the documenation on ipsget()

(See /pkg/rsi/local/libao/phil/ips/a_ipsexample.pro)


IPSCATINP - INPUT AN IPS CATALOG FILE

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipscatinp - input an ips catalog file
SYNTAX: n=ipscatinp(jd,srcI,fname=fname)
ARGS:
jd    : double julian date to use for eps,rsun,hlat,rise,set computations.
KEYWORDS:
fname   : string    name of file to input. def: 'idl/data/ips_strong_to_medium_source.dat'
RETURNS:
 n    : > 0  number of entries in catalog
        0    no data in file
        -1   error.. opening file, or incorrect format
srcI[n]: {}  structure containing the catalog info.

DESCRIPTION:
   input an ips catalog and return info for the requested julian date.
The ips input catalog should have the format:
#J2000     RA(2000)  DEC(2000)   FLUX   SFLUX
 0022+002  002224.7  +001448.6   2.90   1.86
 0032+198  003237.0  +195352.6   4.00   1.83

 The srcname should be the J2000 name without the leading Jj.
 The fluxes are the 327 source and scintillating flux values for the source.

	The routine will compute:
 - lstRise,lstSet: local sideral time rise,set in hours
 - astRise,astSet: ast rise,set times for the source on the specified day.
 - eps : the sun,earth, source angle (in degrees)
 - rsun: the perpendicular distance from sun to earthSource vector (in sun radii).
 - hlat: the heliocentric latitude of the rsun vector.
 
 The stucture contains:
help,srcI,/st                                                                                                         
** Structure <23a2888>, 16 tags, length=160, data length=160, refs=1:                                                       
 BNAME     STRING  'B0019-000'   .. source name (B1950)
 RAB_D     DOUBLE  4.9621558     .. raCoord in degrees (B1950)
 DECB_D    DOUBLE  -0.030337986  .. decCoord in degrees (B1950)
 HHMMSS    STRING  '001950.9'    string version of ra in hhmmss.s
 DDMMSS    STRING  '-000149.2'   string version of dec in ddmmss.s
 FLUX327   DOUBLE  2.9000000     327 flux of source from input cat.
 SFLUX327  DOUBLE  2.9000000     327 scintillating flux
 JD        DOUBLE  2458372.2     Julian date used for computations.
 DATEAST   STRING  'Mon Sep 10 12:00:00 2018' string version of JD (AST)
 RSUN      DOUBLE  65.896268     perpendicular distance sun,earthSrc vector in solar Radii
 EPS       DOUBLE  162.28188     sun,earth,source angle in degrees.
 HLAT      DOUBLE  6.3314646     heliocentric latitude of rsun vector.
 LSTRISE   DOUBLE  23.848367     local siderial rise time (hours)
 LSTSET    DOUBLE  0.93062643    local siderial set  time (hours)
 ASTRISE   DOUBLE  0.95381145    AST rise time for date
 ASTSET    DOUBLE  2.0986467     AST set time for date

Notes:
 The routine limits sources to the AO dec range (-1,38) deg

(See /pkg/rsi/local/libao/phil/ips/ipscatinp.pro)


IPSCATLIST - LIST AN IPS CATALOG

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipscatlist - list an ips catalog
SYNTAX: n=ipscatlist(yymmdd,hhmmss,catlist,duration=duration,jd=jd,catname=catname,hdrLine=hdrLine,$
                      nlinesForHdr=nlinesForHdr,srcI=srcI)
ARGS:
yymmdd: long   date  AST for listing
hhmmss: long   time  AST for listing
KEYWORDS:
JD    : double if supplied then ignore yymmdd,hhmmss, use jd as time for listing
catname: string name of ips catalog to use. def:ips_strong_to_medium_source.dat
nlinesForHdr: long   number of entries before reoutputing the header line.
                     def=40. If <=0 then no header line will be output.
duration: float      duration in hours for listing starting at hhmmss
fout    : string     name for output file. 
RETURNS:
 n    : > 0  number of entries in catalog
        0    no data in file
        -1   error.. opening file, or incorrect format
catlist[n]:strarr   hold formatted catlog
hdrLine: string     header line (if keyword supplied)
srcI[n]: {}         array of structs returned by ipscatinp() (if keyword supplied).

DESCRIPTION:
	input an ips catalog and output a formatted version to stdout. 
The input format should  be:
 srcNm   hhmmss.s ddmmss.s   flux327 scintillatingFlux@327 
a # in column 1 is a comment.

 The input coordinates are assumed  to be J2000. The source Name is the J2000 name without the Jj.
	The output can be used as a catalog file input to cima. The coordinates are translated to
B1950 (including the source name).
 
 The output listing can be used as a catalog input file to cima. The listing will contain:

#srcNm    raB         decB                 date    astRis astSet  eps   rsun  hlat  flux   sflux
B0019-000 001950.9 -000149.2 b 0 Topo v # 20180922 001359 012240 173.3  25.2  17.4   2.9   2.9
B0029+196 002960.0  193720.0 b 0 Topo v # 20180922 233338 022325 157.7  81.8 -42.9   4.0   4.0

 Where:
  astRis,astSet: is the AST rise,set time for the source on this day
  eps          : the sun, earch, source angle (in degrees)
  rsun         : the perpendicular distance sun to earth,source vector (in sun radii)
  hlat         : the heliocentric latitude for the rsun vector.
  flux         : the 327 flux from the input file
  sflux        : the 327 scintillating flux from then input file.

Notes:
 - The program limits sources to the ao dec range: (-1 to +38) degrees.

(See /pkg/rsi/local/libao/phil/ips/ipscatlist.pro)


IPSCMPMINDEX -COMPUTE IPS MODULATION INDEX

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipscmpmindex -compute ips modulation index

SYNTAX: ipscmpmindex,ipsI,ipsspcI,mindexI ,skiptotpwr=skiptotpwr

ARGS:
     ipsI      : struct   from ipsget()
ipsspcI        : struct   from ipscmpspectra()
KEYWORDS:
 skiptotpwr:  if set then skip computing mindex from total power
              pretty worthlless because of rfi.. unless you want to
              compare with spc value as a measure of rfi.. this
              should be the default..

RETURNS:
 mindexI:struct   holds computed modulation index info

DESCRIPTION:

   Compute the modulation index given an ipsI and ipsspcI structs.
The index is computed by subband using the average and then the  median. It 
is also computed for the entire dataset.
  The data is returned in mindexI. 
 When computing the index using the spectra:
   - the median white noise above 65 hz is first removed.
   

IDL> help,mindexI,/st
 SRC      STRING    '1007+075'  ; source name
 YYMMDD   LONG      20180817    ; date from on source filename
 HR       DOUBLE    12.441111   ; hour of day (ast) when observed
 NSB      LONG      5           ; # of subands with data
 RFCFR    DOUBLE    327.00000   ; rf center frequency
 CFRSB    DOUBLE    Array[nsb]  ; rf center freq each subband of data.
 AVGALL   FLOAT     0.454571    ; mindex all data averaging spectra
 MEDALL   FLOAT     0.402009    ; mindex all data using median of spectra
 SPCAVG   FLOAT     Array[nsb, 2] ; mindex by subband and pol using average
 SPCMED   FLOAT     Array[nsb, 2] ; mindex by subband and pol using median
 TP       FLOAT     Array[nsb, 2] ; mindex by suband and pol using 
                                    total power time series.
SEE ALSO:
 ipsget,ipscmpspectra
history
 16apr19:
  - when computing avg all i was using polA twice. fixed to use polA and B
  - when computing white noise to subtract (above 65hz), use the median rather than the mean
  - when integrating spectra, don't include dc, go out to 2 channel before 60hz...

    so radar power doesn't corrupt it so much

(See /pkg/rsi/local/libao/phil/ips/ipscmpmindex.pro)


IPSCMPSPECTRA -COMPUTE IPS SPECTRA

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipscmpspectra -compute ips spectra

SYNTAX: ipscmpspectra,ipsI,reqFreqResHz,ipsSpcI,noblkoff=noblkoff

ARGS:
     ipsI    : struct   from ipsget()
 reqFreqResHz: float    requested freq resolution of the output spectra
KEYWORDS:
noblkoff:    int        added for comet (no equal off). The comet off processing
                        will dump the median off value in a modified ipsI.. 

RETURNS:
 ipsspcI:  struct  holding spectral info

DESCRIPTION:

   Compute the spectra for the total power time series. The user passes in
the ipsI struct returned by ipsget(). The user also supplies the freq resolution
desired for the output spectra. This will determine how many points are used
for each spectra, and how many spectra are returned.
 The following spectral computations are done.
 - for the description below, assume .1 hz resolution, 1 ms sampling, an 80Mhz band
   and 10Mhz subbands.
 - there are 180,000 time samples in 180 secs
 - there are 18 10 second blocks of data
 - there are 10000 samples in a  10 second block
 - for each of the 8 10 Mhz subbands:
    - compute the 18 x 5000 pnt spectra for polA and polB (The input is real so the other 
      5000 points are symmetric.
    - then compute the average of the 18 spectra as well as the median of the 18 spectra.

The data is returned in ipsspcI:
IDL> help,ipsspcIlbwar,/st
 SRC      STRING  'B1005+077'        .. source name
 YYMMDD   LONG    20180824           .. date when on started                                                 
 HR       DOUBLE  12.239444          .. hour of day (ast) when on started                                              
 SPCLEN   LONG    5000               .. length of spectra                                             
 SPCFREQ  DOUBLE  Array[5000]        .. freq array for the spectra                                                     
 NSPC     LONG    18                 .. number of spectra (# of 10 sec blocks)
 NSB      LONG    8                  .. number of frequency sub bands                                          
 BWSB     FLOAT   10.0000            .. bandwidth MHz each sub band                                               
 RFCFR    DOUBLE  1160.0000          .. rf center freq MHz enter band                                              
 RFCFRSB  DOUBLE  Array[8]           .. rf center freq MHz each sub band
 SPC      FLOAT   Array[5000, 18, 2, 8] .. the spectra [spcLen,# 10secBlks,Npol,NfreqSubBands]
 SPCAVG   FLOAT   Array[5000, 2, 8]  .. average of 18 spectra each pol, subband
 SPCMED   FLOAT   Array[5000, 2, 8]  .. median of  18 spectra each pol, subband                                

SEE ALSO:
 ipsget
 27aug18 normalize to on-off
 17may19 change normalization to tpoff for each band
 18may19 switch normalization to area under the spectrum 0.. 55 hz.. 
             ( to stay away from 60hz and radars..

(See /pkg/rsi/local/libao/phil/ips/ipscmpspectra.pro)


IPSDBINP - INPUT THE SAVE FILE HOLDING THE IPS DB INFO

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsdbinp - input the save file holding the ips db info
SYNTAX: n=ipsdbinp(ipsDbAll)
ARGS:
RETURNS:
n         : int  number of entries i ipsDball
                 -1 error.. couldn't find save file
ipsDball[n]: {}  array of structs holding the measured ips info
                 ipsdball is also loaded into the cmips common block by this routine.
DESCRIPTION:
   Input the measured ips info for the measured sources.
After each ips observation, a set of idl.sav files (one for each receiver)  are generated
containing the processed data: total power, on,off diff's, spectra, etc..
	Each save file can be up to 5-10 GBytes. Accessing all of this data is pretty slow.
 	The routine ipsloadallsav() is used to scan all of the save files and extract a
subset of the data into ipsdball[n]. This subset is stored in a single save file (ipsdb.sav).
	ipsdbinp() will input the save file which the user can then query.

	The ipsdball[] struct contains 1 entry for each freq block of a 3min on,off pattern
There can be up to 7 freqblocks for each pattern.
The frequency blocks per pattern  by receiver are:
 327  1  block  with 5 10 Mhz measurements
 lbw  7  blocks each block has 8 10 Mhz measurement
 sbw  7  blocks each block has 8 10Mhz measurements
 sbh  7  blocks each block has 8 10Mhz measurements
 cband 7 blocks each block has 8 20Mhz measuremnets
 xband 7 blocks each block has 8 20Mhz measuremnets
 The ipsdball  stucture contains:
IDL> help,dbsrc,/st
 SRC        STRING  'B0518+165'  src name
 PROJID     STRING  'a3285'      project id for mock datafile
 YYMMDD     LONG     20190526    yyyymmdd date
 HR         FLOAT      12.5833   hour of day (ast)
 AZ         FLOAT      274.703   azimuth deg
 ZA         FLOAT      14.0337   zenith angle deg
 PATID      LONG     914600021   patternId for this observation
 RCVNUM     INT              5   reciever number
 RFCFR      FLOAT      1160.00   rf center freq for this block [MHz]
 NSUBBAND   INT              8   number of subbands used
 SUBBANDCFR DOUBLE    Array[8]   center freq of each subband
 EPS        DOUBLE   17.014711   epsilon deg (earch,sun,src angle)
 HLAT       DOUBLE  -21.648333   helio latitude
 RSUN       DOUBLE   63.743139   distance sun to Perp to observation (solar Radii)
 SRCDEFL    DOUBLE    Array[2]   source deflection polA,B. units=Tsys
 MINDEXAVG  DOUBLE Array[8, 2]   scintillation index [subbands,pols]. average the 10sec spectra
 MINDEXMED  DOUBLE Array[8, 2]   scintillation index [subbands,pols]. median of 10sec spectra
 NMJPEG     STRING 's049_20190526_B0518+165_1160.jpg' Name of jpeg file (used to make movie)
 FON        STRING 'a3285.20190526.b0s1g0.02000.fits' name mock file for first board
 SAVFILENM  STRING 'ips_190526lbw.sav'                name of savefile with raw data.

EXAMPLE:
idl
 @phil
 @ipsinit
istat=ipsdbinp(ipsdball)
srcNm='B0518+165'
n=ipsdbsrc(srcnm,ipsAlldb,dbsrc,/verb)

NOTES:
 - the current data set comes from mock a3285 file observations.

(See /pkg/rsi/local/libao/phil/ips/ipsdbinp.pro)


IPSDBSRC - EXTRACT IPS INFO FROM IPS DATABASE

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsdbsrc - extract ips info from ips database
SYNTAX: n=ipsdbsrc(srcnm,dbsrc,rcvnum=rcvnum,verb=verb)
ARGS:
srcnm   : string  source name. (these are normally B1950 names) of sources
                  we have observed
KEYWORDS:
rcvnum  : int     receiver number: 1-327,5=lbw,7=sbw,8=sbh,9=cband,11=xb
                  If supplied then limit to srcname and this receiver
verb    : int     if set then list the entries in the returned array.
RETURNS:
 n      : int      number of entries found
dbsrc[n] : {}     stuct array holding the observed ips info for the src.

DESCRIPTION:
   Scan the ips database and return a subset holding the requested source.
If rcvnum is also supplied, it will limit this to src and receiver. The
data is returned in dbsrc[] with n being the number of entries.
	The database info is stored in the cmips common block. The common block
is loaded by ipsdbinp(). If the common block has not been loaded, this routine
will call ipsdbinp().

	An entry will be a single frequency block. For each receiver the number of
frequency blocks for a single observation are:
 327  1  block has 5 10 Mhz measurements
 lbw  7  blocks. each block has 8 10 Mhz measurement
 sbw  7  blocks. each block as 8 10Mhz measurements
 sbh  7  blocks each block has 8 10Mhz measurements
 cband 7 blocks each block has 8 20Mhz measuremnets
 xband 7 blocks each block has 8 20Mhz measuremnets
 The dbSrc stucture contains:
IDL> help,dbsrc,/st
 SRC        STRING  'B0518+165'  src name
 PROJID     STRING  'a3285'      project id for mock datafile
 YYMMDD     LONG     20190526    yyyymmdd date
 HR         FLOAT      12.5833   hour of day (ast)
 AZ         FLOAT      274.703   azimuth deg
 ZA         FLOAT      14.0337   zenith angle deg
 PATID      LONG     914600021   patternId for this observation
 RCVNUM     INT              5   reciever number
 RFCFR      FLOAT      1160.00   rf center freq for this block [MHz]
 NSUBBAND   INT              8   number of subbands used
 SUBBANDCFR DOUBLE    Array[8]   center freq of each subband
 EPS        DOUBLE   17.014711   epsilon deg (earch,sun,src angle)
 HLAT       DOUBLE  -21.648333   helio latitude
 RSUN       DOUBLE   63.743139   distance sun to Perp to observation (solar Radii)
 SRCDEFL    DOUBLE    Array[2]   source deflection polA,B. units=Tsys. This is the median
                                 value for the nsubband bands.
 MINDEXAVG  DOUBLE Array[8, 2]   scintillation index [subbands,pols]. average the 10sec spectra
 MINDEXMED  DOUBLE Array[8, 2]   scintillation index [subbands,pols]. median of 10sec spectra
 NMJPEG     STRING 's049_20190526_B0518+165_1160.jpg' Name of jpeg file (used to make movie)
 FON        STRING 'a3285.20190526.b0s1g0.02000.fits' name mock file for first board
 SAVFILENM  STRING 'ips_190526lbw.sav'                name of savefile with raw data.

 When processing a 3 minute observation, we compute 18 10 second spectra.
The mindexxxx (it's a scintillation index, not modulation index:) is computed from the
average of these spectra as well as from the median of these spectra (to check for rfi).

 The /verb option will output a line for each entry kept:
 n=ipsdbsrc('B0518+165',dbsrc,/verb,rcvnum=5)                                                  
yyyymmdd freq    srcDfl    mindexa,mindexb eps   hlat                                              
20190526 1400  2.31  2.21   0.033   0.033  17.0 -21.6                                              
20190604 1400  2.08  2.10   0.067   0.066   9.6 -41.7                                              
 ...
 srcDfl is the source deflection of pola,polB in units of Tsys
 mindexa,mindexb - are the scintillation indices for polA,polB
                   I've takent he median of the nsubband values. 

(See /pkg/rsi/local/libao/phil/ips/ipsdbsrc.pro)


IPSGET - INPUT IPS ON,OFF PAIR, REMOVE RFI AND COMPUTE TOTPWR

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsget - input ips on,off pair, remove rfi and compute totpwr

SYNTAX: istat=ipsget(fon,foff,ipsI,fnmI=fnmI,bon=bon,boff=boff,bfitOff=bfitOff,harmOrder=harmOrder,$
                     bwsubband=bwsubband,spikeClipLev=spikeClipLev,verb=verb)

ARGS:
     fon: string  filename for on source (or use fnmI[2]
    foff: string  filename for off source(or use fnmI[2]
KEYWORDS:
     fmI[2]  : struct   output from masfilelist(). You can use this instead of
                        fon,foff to specify the 2 files
  harmOrder  : long     harmonic order to use when fitting the median off bandpass.
                        def=11
 bwSubBand   : float    bandwidth for subband. For each mock band the spectra
                        will be split into totBw/bwSubBand bands and the total power
                        will be computed separately for each subband. def=10Mhz.
  spikeclipLev: float   When removing narrow spikes from the total power, any spikes
                        greater than sig*spikeClipLev will be interpolated across.
                        . <= 0 --> no spike clipping. def: 4. sigma
     verb: int     if set then print out progress of routine

RETURNS:
 ipsI[]: struct   holding processed info
 istat: 1 ok
      :-1 i/o error..bad hdr,files, not found, etc..
 bon[]: struct    if supplied, the raw spectra from the on  will be returned here
 boff[]:struct    if suppied the raw spectra from the off are retured here
 bfitOff:struct   if suppied then the fit the the median off spectra. This is used
                  for the bandpass correction.


DESCRIPTION:

   Ips datataking uses a 3 min, 3 min off position switch pattern sampled at 2millisecs. This routine
will read in the on and off scans recorded by the mock spectrometers. After input,
 a bandpass correcttion is done, rfi by frequency channel is removed, total power is computed for
 10 Mhz subbands of each mock band, and any short durration spikes are removed.

	Typical sampling is 2 ms for 327,sbw,sbh, and cband. 1 millisecond sampling is used for lbw.
You should use a computer with lots of memory (> 16Gytes) since the on and off file can be
 700Mbytes (1400 MBytes for lbw).
  After input, the processing steps are:
1. compute the median spectrum for the off position.
2. do a robust harmonic fit to the off median spectra. Use the fit for the 
   bandpass correction. The default is an 11 order harmonic fit.
3. divide each of the on,off spectra by the bandpass correction.
4. remove rfi from the spectra. do this for  polA and polB of the on and off spectra:
   - for the description below assume 7 80Mhz bands, 2048 freq channels, 180sec  on,offs and 1ms sampling.
   - process a 10 second block of data at a time (do this 18 times to get 180 secs)
      - compute rms/mean by freq channel. This gives an 2048 array of rms's the freq channels.
      - do a robust linear fit to this rms (throw out 3 sigma outliers, and refit).
      - mark freq channels more than 3 sigma from the fit as not useable.
      - Now split the 80Mhz band into 8 10 Mhz bands. 
      - for each 10 Mhz band (256 freq channels since 2048/8 =256):
         -compute the total power for each 2ms sample. exclude points marked by rms fit as unusable.
          normalize by the number of channels summed.
          The bandpass correction ensures that the included points are weighted correctly.
   - repeat the above for the 18 10 second blocks.
     We end up with total power  arrays of:  blkIOn[18].tp[10000Samples,2pols,8 10mhz bands]
                                       and: blkIOff[18].tp[10000Samples,2pols,8 10mhz bands]
   - Then remove narrow time spikes (1 time sample) from the total power arrays.
     - for each 10Mhz, 180sec tp array subtrace off a median filtered copy (length 3 samples)
     - convolve this with a  [-.5,1.,-.5] kernel (to accentuate spikes of length 1).
     - remove any spike above 3 sigma by interpolating across it. This only applies to narrow spikes.
   - then compute the source deflection: tpDif[90000,2,8]=tpOn - tpOff
   - return the info in the ipsI[]  (it is an array of structures)
The data returned includes:
IDL> help,ipsI,/st
** Structure <1cbacf8>, 21 tags, length=35759736, data length=35759719, refs=1:                      
   SRC          STRING  'B1005+077'   .. used                                                            
   FON          STRING  '/share/pdata1/pdev/a3285.20180824.b0s1g0.01200.fits'                   
   FOFF         STRING  '/share/pdata1/pdev/a3285.20180824.b0s1g0.01300.fits'                   
   NREC         LONG    180           .. number of records/seconds in on or off                                         
   SMPREC       LONG    1000          .. 1 second records written to disc.                                             
   SMPTOT       LONG    180000        .. number of 1 millisecond samples                                                
   BWTOT        DOUBLE  80.000000     .. total bandwidth 1 mock box                                                    
   BWSUBBAND    FLOAT   10.0000       .. subband bandwidth (MHz)                                                    
   NSUBBAND     LONG    8             .. number of sub bands                                               
   RFCFR        DOUBLE  1160.0000     .. rf center frequency of band                                                   
   RFFREQ       DOUBLE  Array[2048]   .. rf frequency for the input spectra                                             
   CFRSUBBAND   DOUBLE  Array[8]      .. rf center frequency each subband
   HDR          STRUCT  -> MASFHDR Array[1] .. the header of the 1st record of the on fits file
   BLKLEN       INT     10            .. data split into 10 sec blocks for rfi removal
   NBLKS        LONG    18            .. number of blocks in 180 secs          
   BLKION       STRUCT  ->  Array[18] ..total power info processing on source
   BLKIOFF      STRUCT  ->  Array[18] ..total power info processing off source
   SPIKECLIPLEV FLOAT   4.            .. nsigma for total power narrow  spike clipping
   NCLIPON      LONG    Array[2, 8]   .. # samples clipped each subband [2pol,8subbands] on src
   NCLIPOFF     LONG    Array[2, 8]   .. # samples clipped each subband off src
   TPDIF        FLOAT   Array[180000, 2, 8] ..Total power after removing narrow spikes

 The BLkIon,off structures contain:
IDL> help,ipsI.blkIon,/st                                                                            
** Structure <1c21a58>, 6 tags, length=672836, data length=672836, refs=3:                           
 I0REC     INT       0                 record # start 10 sec block (cnt from 0)
 I1REC     INT       9                 record # end 10 sec block                      
 RMSBYCHAN FLOAT     Array[2048, 2]    rms/mean computed for each freq chan polA,B
 MASK      FLOAT     Array[2048, 2]    good freq channels.1=good, 0- skip
 NGOOD     LONG      Array[2, 8]       number of good channels 2pol, 8 subbands
 TP        FLOAT     Array[10000,2,8]  total power 10000 samples, 2pol, 8 subbands
                                       computed after exclude bad frequency channels.                     

 The raw spectral data read from the fits file is optionally retured in the bon and/or boff keywords
 (if supplied by the user). This is the standard mas routine format (see masdoc).
 
 Setting the verb keyword to a non zero value will cause the routine to output its progress while it is running.

EXAMPLE:
 	To use these routines you should:
 idl71
 @phil
 @ipsinit
 
 The following example would process data:
  -  taken by project a3285  on 14aug18
  - using a single mock board at 327 mhz..
  - the idl code would be:
  proj='a3285'
  yymmdd=20180814
  n=masfilelist(yymmdd,fnmI,proj=proj,/appbm,grp=0,band=1) .. scan for all files for this dataset
;  .. if you didn't include bm0 in the data acquisition, you might have to include
;  .. the bm=  parameter (bm is same as box in the cima menu).
  n=masfilesum(flist,sumI,fnmI=fnmI,/list)                 .. read header info for each file
  ; loop printing info on each file to see what's there
  for i=0,n-1 do print,i,sumI[i].num,sumI[i].h.object,sumI[i].h.obsmode,sumI[i].h.scantype,$
                        format='(i2," fileNum:",i4," src:",a10," patType:",a," scanType:",a)'
0 fileNum:   0 src: B0758+143 patType:ONOFF scanType:ON
 1 fileNum: 100 src: B0758+143 patType:ONOFF scanType:OFF
 2 fileNum: 200 src:  0952+284 patType:ONOFF scanType:ON
 3 fileNum: 300 src:  0952+284 patType:ONOFF scanType:OFF
 4 fileNum: 400 src:     3C237 patType:ONOFF scanType:ON
 5 fileNum: 500 src:     3C237 patType:ONOFF scanType:OFF
 6 fileNum: 600 src:  1011+064 patType:ONOFF scanType:ON
 7 fileNum: 700 src:  1011+064 patType:ONOFF scanType:OFF
 8 fileNum: 800 src:  0936+043 patType:ONOFF scanType:ON
 9 fileNum: 900 src:  0936+043 patType:ONOFF scanType:OFF
10 fileNum:1000 src:  0947+000 patType:ONOFF scanType:ON
11 fileNum:1100 src:  0947+000 patType:ONOFF scanType:OFF

;.. nlet's look at file num 600,700 (indices 6,7)

 fnmIL=fnmI[6:7]   ; uses indices 6,7
; now call the routine.. include bon,boff keywords so the raw data is also returned.
 istat=ipsget(fon,foff,ipsI,fnmI=fnmIL,bon=bon,boff=boff,/verb)
; plot the total power  for on and off (polA and B still separate) for 1 subband
!p.multi=[0,1,2]
 isb=0   .. the first 10 Mhz subband
 ver,800,2000
; polA
 ipol=0
 hor
 plot,ipsI.blkIOn.tp[*,ipol,isb],tit='toal power. polA'
 oplot,ipsI.blkIOff.tp[*,ipol,isb],col=colph[2]
; polB
 ipol=1
 hor
 plot,ipsI.blkIOn.tp[*,ipol,isb],tit='toal power. polB'
 oplot,ipsI.blkIOff.tp[*,ipol,isb],col=colph[2]

; You could also plot the source deflection for the 1st subband
 ipol=0
 plot,ipsI.tpDif[*,ipol,isb],tit='Source deflection. white polA, red polB'
 ipol=1
 oplot,ipsI.tpDif[*,ipol,isb],col=colph[2]
SEE ALSO:
 ipsremovespikes, mas routines

(See /pkg/rsi/local/libao/phil/ips/ipsget.pro)


IPSPOSITION - COMPUTE SOURCE,SUN ANGLES FOR IPS OBSERVATION.

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsposition - compute source,sun angles for ips observation.

SYNTAX: ipsposition,srcNm,jd,ra,dec,posI,b1950=b1950
                     ipsI=ipsI

ARGS:
   srcNm: string  name of source.
     jd : double  julian date (can be a scalar or vector)
                  if a vector then the ra,dec will be evaluated at
                  multiple dates.
    ra  : double  ra position in degrees (scalar) def to J2000
   dec  ; double  dec position in degrees (scalar)def to J2000
KEYWORDS:
 b1950 :      int      if set that ra,dec is in B1950 coord. Default is J2000
 ipsI  : {}       ipsI struct returned by ipsget(). If supplied,
                  then take all the needed info from here.. Ignore the 
                  the other args.

RETURNS:
 posI[]: struct   holding processed info

DESCRIPTION:
   For a given ra,dec, compute the position info for the specified julian dates.
The user specifies:
  - the source name
  - ra (in degrees), deg (in degrees). the default coord sys is J2000
       if the coord are B1950 then you need to set the keyword /b1950
	The routine computes position info for the source coord at each of the
specified julian dates.

	A second way to call the routine is to use the ipsI=ipsI   keyword.
In this case, you pass in an array of ipsI structures (from ipsget(), or from the
 .sav files). The routine will then compute the position info for each entry
in the ipsI array. The other input arguements are ignored (although dummy values need to be 
supplied).

The position Info includes:
   - The source position as:
       - input ra,dec
       - ra,dec of date
       - geocentric eclipitc position of date
   - earth sun 3d vector in geocentric ecliptic coord of date
   - sun,earth,src angle (epsilon)
   - radial distance sun to perpendicular to earth,source vector
   - the heliocentric latitude of above radial distance.
 help,posI,/st
  _D --> degrees
   C --> current coordinate
   E --> geocentric ecliptic system
  ** Structure <242ba38>, 15 tags, length=152, data length=146, refs=1:
 SRC         STRING 'B1005+077'  .. source name
 JD          DOUBLE  2458359.1   .. julian date
 DATEAST     STRING 'Tue Aug 28 11:29:04 2018' ast time start obs.
 SRCRA_D     DOUBLE  170.11580  source ra degrees input. Seee b1950 for j or B
 SRCDEC_D    DOUBLE  14.348508  source deg degrees input
 B1950       INT     0          if  0 the srcra,dec is J2000. if 1 then B1950
 SRCRAC_D    DOUBLE  152.24507  source ra Current coordiantes. deg
 SRCDECC_D   DOUBLE  7.4127896  source dec Current coordinates. deg
 SRCLONGEC_D DOUBLE  151.57365  source longitude geocentric ecliptic coord
 SRCLATEC_D  DOUBLE  -3.7453532 source latitude geocentric eclipitic coord
 SUNV_EC_AU  DOUBLE  Array[3]   Earth to sun vector AU, geocentric eclipitic coord
 SUNDISTAU   DOUBLE  1.0101170  sun earth distanance in AU
 EPS_D       DOUBLE  5.2639184  epsilon.sun,earth, source angle
 HLAT_D      DOUBLE  -47.796915 Helicentric ecliptic latitude of rsun deg.
 RSUN        DOUBLE  19.927393  distance sun to perpendiculr to earch source vector
                                in sun radii
 The routine uses the idl goddard routines: sunpos(),xyz(), and precess

EXAMPLES:
  astToUtc=4./24d        ; ao is 4 hours behind utc
; pass in parameters with 2 julian days
  jd=[julday(8,28,2018d,12,0,0),julday(8,29,2018d,12,0,0)] + 4d/24
  src='B1005+077'
  ra =170.11580
  dec=14.348508
  ipsposition,src,ra,dec,posI
  help,posI
  POSI            STRUCT    = ->  Array[2]
; get position info for all of the  sbh patterns on 28aug18

  restore,'/share/ips/a3285/ips_180828sbh.sav',/verb
% RESTORE: Portable (XDR) SAVE/RESTORE file.
% RESTORE: Save file written by phil@gpuserv0, Fri Aug 31 16:49:35 2018.
% RESTORE: IDL version 7.1 (linux, x86_64).
% RESTORE: Restored variable: YYMMDD.
% RESTORE: Restored variable: NPAT.
% RESTORE: Restored variable: NBAND.
% RESTORE: Restored variable: IPSISBHAR.
% RESTORE: Restored variable: IPSSPCISBHAR.
 help,ipsIsbhAr
 IPSISBHAR       STRUCT    = ->  Array[42] ; 42 patterns

	ipsposition,'',0,0,0,posI,ipsI=ipsIsbhar
  help,posI
  POSI            STRUCT    = ->  Array[42]
Notes:
   - Ra's are expected in degrees, not hours
   - The ecliptic value are geocentric (except for heliocentric latitude). 
   - The first calling method (without ipsI= keyword) takes 1 ra,dec and 1 or more
     julian dates jd[]. It computes the position info for the single ra,dec at
     multiple times
   - The second calling method (using ipsI=ipsI) takes everything from the ipsI struct
     It uses the ra,dec, jd from each entry for the computation. So you can have
     multiple ra,decs and jds..

(See /pkg/rsi/local/libao/phil/ips/ipsposition.pro)


IPSREMOVESPIKES -REMOVE NARROW SPIKES FROM TOTAL POWER

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsremovespikes -remove narrow spikes from total power

SYNTAX: nremoved=ipsremovespikes(tpIn,tpOut,indspikes=indspikes,nsig=nsig,medlen=medlen)

ARGS:
     tpIn[n]: float      total power data
    nsig    : float      to clip spikes. def: 4 sigma.
   medlen   : long       length of median filter to use for subraction
KEYWORDS:

RETURNS:
 nremoved: long          number of spikes removed and interpolated over
 tpOut[n]: float         data with spikes removed
 indspikes[nremoved]: long   indices for spikes that were removed

DESCRIPTION:

   Try to remove narrow spikes (1 channel wide) in total power data. The
program will:
1. create a median filtered verision of the input, then subtract from tpIn
    tpDif=tpIn -median(tpIn,medLen) ..default medlen =3
3. convolve a kernel of [-.5,1.,-.5] with tpDif. this will enhance narrow spikes
4. flag all points greater than the clipping level (def=.04). The on and off
   have been scaled to median(boff) so the unit values are around 1.
5. for each flagged point, interpolate a new value using the adjacent points
   (this does not check to see if an adjacent point is also flagged ...)
6  return tpOut with the spikes removed and well as the number of spiked removed.
  . optionally return the  indices for the spikes in indspike[]
SEE ALSO:
 ipsget

(See /pkg/rsi/local/libao/phil/ips/ipsremovespikes.pro)


IPSSCANPATTERNS - SCAN FOR IPS PATTERNS FOR A DAY

[Previous Routine] [Next Routine] [List of Routines]
NAME:
ipsscanpatterns - scan for ips patterns for a day

SYNTAX: n=ipsscanpatterns(projId,yyyymmdd,patI,verb=verb)

ARGS:
 projId      : string  project id to scan for
 yyyymmdd    : long    date to scan
KEYWORDS:
 verb        : int     if set then print what we found to stdout

RETURNS:
 istat:  int :  < 0 error
                >=0  number of patterns we found.
 patI[n]     : {}    array of structs holding info on each pattern.

DESCRIPTION:

   Scan for on,off patterns taken by projId on the date yyyymmdd. Return
the info in patI

The data is returned in patI is
IDL> help,patI,/st
 SRCNM    STRING    'J0936+043'  .. source name
 BM       INT              0     .. mock beam/box 0..6
 GRP      INT              0     .. mock grp .. 0 or 1
 NSCANS   INT              2     .. # of scan=2  on,off
 PATID    LONG         823600018 .. pattern id number
 FILELIST STRING    Array[2]     .. on,off filenames
 SUMI     STRUCT    ->  Array[2] .. hold summerinfo from on and off scan

 IDL>  help,patI.sumI,/st
** Structure , 13 tags, length=1256, data length=1241, refs=2:
 OK       INT              1
 FNAME    STRING    '/share/pdata1/pdev/a3285.20180824.b0s1g0.00200.fits' filename
 FSIZE    LONG64                 744223680 .. bytes in  file
 NROWS    LONG               181  number of rows/recs (last is a partial rec)
 DUMPROW  LONG               500  spectra per row
 DATE     LONG          20180824  .. date  the next 5 are parsed from the file name
 BM       INT              0      .. beam
 BAND     INT              1      .. band
 GRP      INT              0      .. grp
 NUM      LONG               200  .. file number
 H        STRUCT    -> MASFHDR Array[1] .. fits header
 HMAIN    STRUCT    -> PDEV_HDRPDEV Array[1] .. mock device header
 HSP1     STRUCT    -> PDEV_HDRSP1 Array[1]  .. mock fpga header

 Gotchas: 
 - routine setup to only look at band 1.. all single pixel receivers only
   use band 1.. alfa uses band0 and band1 .. this  will skip 1/2 of alfas files
 - only looks for grp0 files. if both mock groups are used (piggy back observing)
   then it will skip grp 1.

 SEE ALSO: maspatfind

(See /pkg/rsi/local/libao/phil/ips/ipsscanpatterns.pro)


IPSSRCINFO - PRINT AND RETURN MEASURED RESULTS OF IPS DATA

[Previous Routine] [List of Routines]
NAME:
ipssrcinfo - print and return measured results of ips data
SYNTAX: n=ipssrcinfo(src,srcI,yymmdd=yymmdd,rcvnum=rcvnum)
ARGS:
 src:  string  name of source to search for.
KEYWORDS:
rcvnum: int     if supplied the limit returned data to specified
                receiver. The rcvnums are:
                 (1-327),(5-lbw),(7-sbw),(8-sbh),(9=cband),(11-xband)
yymmdd: long    if supplied then return all sources from date rather
                then all dates for src. In this case params
                src is ignored
RETURNS:
N   :  int     number of separate measurements found
srcI[N]: {}    array of struct holding the info
DESCRIPTION:
 	Search the ao idl ips database for all measurements done for
a specified set of parameters. The normal usage is for the 
user to specify a source name. It will then return all measurements
for that source. 
	If the keyword rcvnum= is supplied then measurements will be 
limited to the specified receiver.
 	If the keyword yymmdd= is supplied, then the src parameter
is ignored, and all measurements for the specified data are returned
(the src parameter can be ''). 

Note:
	The info printed/returned is a subset of the recorded data. This
routine is helpful when setting up observations. It will tell you
whether an observation at a particulur epsilon and frequency will give
are reasonable result. For the full data set use ipsdbsr().

	The returned values are:
n - the number of entries found
srcI[n] - struct holding some of the source info. it contains:

help,srcI,/st
Structure <12ed7a8>, 6 tags, length=64, data length=64, refs=1:
SRC      STRING    'B2230+114'       ..the source name
YYMMDD   LONG          20200318      ..the date of the observation:yyyymmdd
RFCFR    FLOAT           1400.00     ..rf center freq for one of the 80 Mhz bands
SRCDEFL  DOUBLE    Array[2]          ..polA,polB  The median srcDeflection/Tsys 
                                       for the nsubands (usually 80Mhz)
EPS      DOUBLE           23.620274  .. earth,sun,source angle for measurement
MINDEX   DOUBLE    Array[2]          .. PolA,polB. median scintillation index for 
                                        the measurement over the nsubands (usually 80 Mhz).
 The data printed to the terminl is:
IDL> n=ipssrcinfo(src,srcI)
freq   eps     srcDefl        mindex   yymmdd
 327  39.81   0.68   0.34 0.057 0.058 20200409
 327  40.65   0.69   0.35 0.055 0.055 20200410
1400  23.62   1.96   1.96 0.019 0.019 20200318
1400  24.19   1.97   1.93 0.017 0.018 20200319
1400  24.79   1.89   2.00 0.015 0.015 20200320

 It is sorted by frequency,, and then epsilon  within frequency

(See /pkg/rsi/local/libao/phil/ips/ipssrcinfo.pro)