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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)