NAME:
atmclp - coded long pulse (ri)
SYNTAX: istat=atmclp(lun,spcBuf1,spcBuf2,nrec=nrec,baudLen=baudLen,
spclen=spclen,txSmpSkip=txSmpSkip,$
dinfo=dinfo,hdr=hdr,dotm=dotm,tmi=tmi
ARGS:
lun: int lun that points at data file
KEYWORDS:
nrec: long number of records to process (default is 1000).
warning.. this is the number of records, not ipps
unless ipps/buf = 1.
baudlen: float By default the baudlen is taken from the header. Some
of the datataking files have set the baudlen
equal to the codelength. If this is found then the
baudlen is set by default to 1 usec. You can force the
the baudlen by setting this keyword (the units are usecs).
spclen : long The length of the spectra to do. By default it is
rounded up to the next power of 2.
txSmpSkip:float The number of usecs to skip before taking the first
tx sample. This takes in account the filter delay for the
tx samples. The default is 2 usecs.
dotm : if set then do detailed tming.. if not, just do
total times.
RETURNS:
spcbuf1[spclen,nhghts]: float the averaged spectra vs heights for the
first fifo (normally ch1 for dual beam).
spcbuf2[spclen,nhghts]: float the averaged spectra vs height if two
fifos were used (normally gr for dual beam).
dinfo : {} Structure holding info that was used in the computation
(see below)
hdr : {} header from the first record averaged.
tmI : {} Timing info (see timing info below for a description). If
dotm=1 then you get detailed info. If not, then just the totals.
DESCRIPTION:
Input and process coded long pulse data taken in raw data mode with
the radar interface. It will read nrec records (default 1000) starting
from the current location on disc (pointed to by lun). It will read the
requested number of records into a buffer (so don't make it too large). It
will stop reading prematurely if it finds a data record of a different
type (mracf, power, tp).
For each ipp the processing step is:
1. extract the tx samples using txSmpSkip to determine which
tx Sample to start on.
2. compute the complex conjugate of these samples.
3. for each height:
- multiply the data by the code.
- zero extend to spclen
- fft the data
- compute power and accumulate
- compute how many samples to skip to get to the next starting
height (usually baudlen/sampLen).
4. After accumulating divide by the number of accumulations;
5. Shift each spectra so that dc is in the center (for a 4096 length
xform, shift right by 2048)
6. fill in the dinfo structure with the info that was used for the
computation. Dinfo contaings:
GWUSEC FLOAT 0.200000
BAUDLENUSEC FLOAT 1.00000
CODELENGW LONG 2500
NHGHTS LONG 601
HGHTSTEPGW LONG 5
IPPAVGED LONG 1000
SPCLEN LONG 4096
TXSMPSKIP LONG 10
NUMFIFO INT 2
7. also return the header from the first ipp used in the variable hdr.
Interesting info is in the sps portion of the header (although the
baudlen may be incorrect).
WARNING:
This routine was tested on one set of data (dual beam, 1 ipp/buf, 5Mhz bw, 1 usec
buad). The results looked pretty good (ie we saw the plasma line..), but a detailed
comparison with the asp version has not been done.
TIMING INFO:
Timing info can be returned if tmi=tmi keyword is provided. The tmI
structure holds times for different parts of the code. Each time
structure has the number of times the code was timed, the total sum
and the total sum of squares. The routine printtmI,tmi can be used
to print this info out.
An example output is:
IDL> printtmi,tmi
tmTot:675.096 read: 3.987
CODECONJ Ntot: 1000 tmTot: 0.4328 avgTm0.000433 sig:0.000006
CODEM1 Ntot:601000 tmTot: 51.1495 avgTm0.000085 sig:0.000005
FFT1 Ntot:601000 tmTot:122.8049 avgTm0.000204 sig:0.000008
ACCUM1 Ntot:601000 tmTot:135.5683 avgTm0.000226 sig:0.000006
CODEM2 Ntot:601000 tmTot: 51.2731 avgTm0.000085 sig:0.000004
FFT2 Ntot:601000 tmTot:122.8776 avgTm0.000204 sig:0.000007
ACCUM2 Ntot:601000 tmTot:135.2315 avgTm0.000225 sig:0.000006
BUF1TOT Ntot:601000 tmTot:322.3068 avgTm0.000536 sig:0.000013
BUF2TOT Ntot:601000 tmTot:322.0795 avgTm0.000536 sig:0.000012
All times are in seconds. the columns are:
totmeasTm total wall time
read: total time for reading data.
sectionNm The section of code that was timed.
Ntot is how many times this section was timed.,
tmTot is the total time for this section.
avgTm is the average time for 1 execution of this code
sig is the sigma for the Ntot measurements
The sections are:
CODECONJ : extract and conjugate the tx samples
CODEM1 : extract hght data, multiply by code, 0 extend.
FFT1 : compute the fft (using fttw)
accum1 : compute the power and accumulate
BUF1TOT : total time for buf1. includes some data manipulation times
not included in the above times
The xxx2 : same as xxx1 but for the second fifo (if present)
totmeasTm : this is just the sum of the measured times. It does not
include the overhead of timing (about 2usecs per call).
The example above had 1000 codeconj so 1000 ipps were processed (10 secs
of data).
You'll notice that computing the power and accumulating is taking
longer than the fft.
For this computer fftw in C took about 44 usecs for a 1k xform. So it
should take about 211 usecs for a 4k. It is averaging 204 so there doesn't
look like there is any speed up time in the fft.
Code speedup would probably come by writing an external module
that was passed a tx rec, hght rec, and then it did the fftw and
power, accumulate, returing the averaged buf. This might give you
50% speed up.
(See /pkg/rsi/local/libao/phil/atm/atmclp.pro)
NAME:
atmdcd - decode atm rawdat data
SYNTAX: nipps=atmdcd(d,code,vdcd,nhghts,dcdh=dcdh,firsttime=firsttime,
use2ndchan=use2ndchan,usecal=usecal,
barkercode=barkercode,codelen88=codelen88,codelen52=codelen52)
ARGS:
d[n]:{rd} array of rawdat structures to decode. These are
returned by atmget().
code[codelen]:float code to use (unless one of the codes has been specified
with one of the code keywords.
KEYWORDS:
dcdH:{} decode structure holding info for future calls
using the same data setup. It speeds up future
calls by caching the transformed code and other info
so they don't need to be compute each time.
firstTime: If set, then this is the first call with the data set.
If provided, it will load dcdh with info to speed up
future calls.
use2ndchan: if set then decode the 2nd chnannel of the data
in d. By default the first channel is always decoded
usecal : if set then include any cal samples in the decoding
barkercode: if set then ignore code[]. Program will generate the
barker code and it use for decoding.
codelen88 : if set then ignore code[]. Program will generate the
88 length code (1 usec meteor observations) and use it
for decoding.
codelen52 : if set then ignore code[]. Program will generate the
52 length code used for the dregion and use it for decoding.
RETURNS:
nipps:long number of ipps decoded
vdcd[nhghts,nipps]:complex the decoded voltages
nhghts :long number of heights decoded..
dcdh:{} structure holding info to use to speed up
calls after the first call. On the first call
this is returned to the user. On succeeding
calls it is just an input parameter.
DESCRIPTION:
Decode a number of ipps from atm rawdat. The user passes in data
from atmget as well as the code to use. The routine will decode the
ipps and pass the decoded voltages back in vdcd.
The program has 3 built in codes. If you set the keyword that corresponds
to one of these codes, then the program will generate the code for you
(and it will ignore whatever you pass in via code). On return, code()
will contain the generated code.
When decoding you need to zero extend and transform the code. If you
are decoding many records with the same kind of data, you only need to compute
the spectrum of the code once. The program lets you do this by:
1. On the first call the user sets firsttime=1 and puts
dcdh=dcdh on the call so atmdcd can return info in it.
2. atmdcd zero extends and transforms the code. It places this and
other info into dcdh which then gets passed back to the user.
3. On subsequent calls the user sets firsttime=0 and passes in dcdh=dcdh.
4. atmget just grabs the cached info from dcdh rather than recomputing
it.
If you change the data type (hghts, codelen,code, txsamples, channel
to decode, then you must set firsttime=1 again to recompute the new
parameters.
Description of dcdh structure:
dcdh.CODELEN LONG 52 Length of code used
dcdh.FFTLEN LONG 512 length of fft used for decoding
dcdh.DCDHGHTS LONG 349 number of fully decoded hghts
dcdh.SMPTX LONG 0 number of tx samples
dcdh.SMPHGHT LONG 400 number of sample heights
dcdh.SMP1IPP LONG 400 number of samples in ipp
dcdh.IPPSREC LONG 4 number of ipps in record
dcdh.INDFIFO INT 1 index into d.(indfifo) for data buf to use
dcdh.FFTSCALE FLOAT 262144. scaling factor for fft (fftlen)^2
dcdh.SPCCODE COMPLEX Array[512] zero exteneded transformed code,conjugated
and then scaled by fftscale.
If this is too confusing, just ignore firsttime= and dcdh=. Everything
will work ok but it'll run a lot slower since it recomputes the code on
each call.
The last codelen-1 heights are not decoded (since all the data is not
present). The returned array holds hghts:
0 through Nhghts-(codelen-1) heights
If a cal is present, then the /usecal keyword will cause the cal
samples to be included in the decoding (this is how the old power profile
program decoded the data).
EXAMPLES:
1. loop calling atmget and atmdcd. Assume that the data was coded
with the 52 length d region code. Don't bother to remove the
digitizer offsets. Assume there are 4 ipps/record and that you
want 256*300 ipps to process:
the data.
firsttime=1
recPerRead=100 ; read 100 records per atmget call
nreads=100 ; read 100 records per atmget call
ip=0L
for i=0,nreads-1 do begin
istat=atmget(lun,d,nrec=recPerRead,/search)
nipps=atmdcd(d,code,vdcd,nhght,/codelen52,firsttime=firsttime,dcdh=dcdh)
firsttime=0 ; after 1st time use dcdh.
if i eq 0 then begin ; first time allocate volt array
ippsTot=nreads*recPerRead*d[0].h.ri.ippsperbuf
vAr=complexarr(nhght,ippsTot)
endif
vAr[*,ip:ip+nipps-1]=vdcd
ipp+=nipps
endfor
The var[ndcdhghts:349,ippsTot:40000.] will use 112 megabytes.
NOTES:
The current scaling for atmdcd() will decode a code of unit height to
a single range of amplitude codelen.
(See /pkg/rsi/local/libao/phil/atm/atmdcd.pro)
NAME:
atmdcdppspcavg - decode, compute pulse to pulse spc, and avg
SYNTAX: nspc=atmdcdppspcavg(lun,secToAvg,spcLen,spcAvg,h=h,
use2ndChn=use2ndChn,useMedian=useMedian,verb=verb,
rectype=rectype,spctoavg=spctoavg)
ARGS:
lun:int lun for file to process. You should have alread opened it.
secToAvg:float seconds of data to average spectra over.
spclen:int length for spectra. This is the number of ipps 1
spectra will cover.
KEYWORDS:
use2ndchan: if set then process the 2nd channel of data (fifo 2). The
default is the first channel.
useMedian: if set then use the median when averaging the spectra.
The default is the mean.
verb: if set then print out when we start each block of spectra.
rectype: string Type of record to process. Use this if different kinds
of records are present in the same file. The rectypes
are defined in atmget(). For dregino profiles use
rpwr88 or r52code. Be careful using this routine with
datasets that are not contiguous in time. The spclen must
divide evenly into the records for a given time block.
spcToavg: int if provided then ignore sectoavg, use this.. works
better when ipp is not a submultiple of 1 sec
RETURNS:
nspc:long The number of spectra that have been averaged.
spcAvg[spcLen,nhghts]:float the averaged spectra. One for each decoded
height.
h:{} The header from the first record input.
DESCRIPTION:
Decode a number of ipps from atm rawdat, compute the pulse to pulse
spectra, and then average these together. The user passes in the lun for
the file to process, the seconds of data to process, and the length of the
transform to use for the spectra.
The routine will process a block of spectra at a time. A block
is defined as the number of spectra that will fit in 700 mb of memory
(or the requested seconds of data if it is smaller).
For each block of spectra:
1. the data is input.
2. the mean value is removed from the complex voltages (dc is removed)
3. the data is decoded
4. the transform along the height direction is done for
tne number of spc in the block
5. power is computeded
6. the spectral power is summed for each height (number of spectra in
the block).
If the usemedian keyword is set then step 6 uses the median value for
each spectral height (over the number of spectra in the block). At the
end the median from each block are averaged (weighted by the number of
spectra in each median).
For a single channel of data (ch or dome), 349 heights, and 1 millisecond
ipps, a block of spectra will be about 2 minutes of data.
EXAMPLES:
average 3 minutes of data:
file='/share/xserve0/aeron5bup/T2166/t2166_21mar2006.023'
openr,lun,file,/get_lun
secToAvg=60*3
spcLen =512
;
; this data set has:
; 1. an ipp of 1040 usec, and 400 samples in the first rcv window.
; 2. the codelen is 52. So there are 400-51=349 hghts after decoding.
; 3. 3 minutes of data is 173077 ipps. This is 339 512 length spectra
; 4. 227 spectra fit in a block of 700Mb . The first block has 227 spectra,
; the second block has 112 spectra
nspc=atmdcdppspcavg(lun,secToAvg,spcLen,spcAvgMed,/verb,/usemed)
help,spcavgmed,/st
SPCAVGMED FLOAT = Array[512, 349]
; make an image of the results
minval=5e-3
maxval=1.6e4
imgdisp,(spcavgMed > minval)
(See /pkg/rsi/local/libao/phil/atm/atmdcdppspcavg.pro)
ATMGET - INPUT 1 OR MORE ATM RECORDS
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
atmget - input 1 or more atm records
SYNTAX: istat=atmget(lun,d,nrec=nrec,rectype=rectype,usercodelen=usercodelen,
raw=raw,scan=scan,sl=sl,search=search,
contiguous=continguous)
ARGS: lun : unit number for file (already opened)
d[] : {atmrec} return data here.
KEYWORDS:
nrec : long number of records to return. The default is 1. If nrec is
greater than 1, the routine will skip over any records that
do not match the requested record type.
Note: The routine will return when: nrec records
of the same type are read, eof is hit, or a record of
the same record type but with a different number of
data elements is found.
rectype : string record type to return. The default is the next
available record. The possible values are:
'pwr' : power profile (ap version)
'clp' : codedlong pulse (ap version)
'mracf': mracf (ap version)
'dspc': dynamic spectra (ap version)
'rawdat': raw dat program (any data program)
'rpwrb': rawDat barker code power profile
'rpwr88': rawDat 88 len code (power profile or d-region use)
'r52code': rawDat 52 length code
'rmracf': rawDat mracf;
'ruser' : rawDat user specified code. You must also
set usrcodelen=usercodenlen keyword to the
length of the code in bauds
'rclp': rawDat coded long pulse
note: rmracf,rpwrb,rpwr88,r52code use the codename to
determine the type of record.
usercodelen: long if rectype is ruser, you must also set usercodelen
to the baud length of the code you want to read
raw: if set then return the data in a single float array.
By default the routine splits the data in each
record into tx,data,noise,cals, etc..
(but rawdat,rpwr,rmracf always return just the data)
scan: long scan to position to before reading the data.
sl[]: {scanlist} array returned from getsl(lun). This can be used
for random access of files (see generic routines: getsl())
search: if set then search for the start of the first header.
use this in if the first rec of the file does not start
with a header.
contiguous: if set then the records returned will be contiguous
in time. This is handy when a single file has different types
of records in it. Suppose the file holds:
200 rclp recs,
1024 'rpwr88' recs
200 rclp recs,
1024 'rpwr88' recs
etc.
If you set rectype='rpwr88' and ask for 2000 records, then
it will return with 1024 records. If contiguous is not set,
then it will read 1024 recs of the first block and then
read the rest from the second 1024 block. This can be a
problem if doing pulse to pulse spectra.
nrecs=1000. If after
RETURNS:
istat : int
1 - requested number of records found
2 - req number of recs not found, hit eof, at least 1 record found
3 - req number of recs not found, hit rec with different data len,
at least 1 record found
0 - hit eof, no data returned
-1 - could not position to scan/rec
-2 - bad header id in hdr
-3 - bad program id in header
-4 - requested rectype='ruser' but forgot to set usercodelen
NOTE:
For now the rmracf,r
DESCRIPTION:
Read the next nrec atm records from lun. Start at the current
record position unless the scan keyword is present. Atm files contain
more than 1 kind of data records (eg. pwr, mracf, clp, etc). By default
the record type returned is determined by the first record of the current
read. The user can use the rectype keyword to specify a particular kind
of record type to return.
The file is left positioned after the last successful record read. If
eof() is hit while trying to complete the current request, the file will
be left positioned at the last successful read of the requested record
type. If the file does not end in this particular record type, then this
will not be at the end of the file.
The data is returned in an array of structures d[nrec]. Each element
of the structure will contain the header followed by the data arrays. By
default the data arrays are split into tx,data,noise,cals, etc. The
raw keyword will return the data in each record as a single float array.
EXAMPLES:
To use this routine:
0. idl
1. make sure the path to aodefdir() is in your path directory:
at ao just enter @phil or
!path='/home/phil/idl/gen' + !path
You can also put this in your IDL_STARTUP file if you like.
2. @atminit .. sets up the paths, headers for atm processing.
3. openr,lun,'/share/aeron2/29jun03qz09.greg',/get_lun
--> get 10 pwr records
istat=atmget(lun,d,rectype='pwr',nrecs=10)
--> get 20 mracf records, return in raw mode without splitting
up the data into tx,data,cals,etc.
istat=atmget(lun,d,rectype='mracf',nrecs=20,/raw)
If you want to see each record of a file try:
rew,lun
scanlist,lun,/verb,/std
rec: 1 scn: 0 grp: 251 Id:clp pos: 0 h,dlen: 444 154112
rec: 2 scn: 0 grp: 501 Id:clp pos: 154556 h,dlen: 444 154112
The position is the byteoffset in the file for the start of the record.
h,dlen are the header,data lengths in bytes.
-------------------------------------------------------------------------
RESTRICTIONS:
1. This routine currently works for the data taken in "raw datataking" mode.
This uses the pc datataking system in raw datataking mode with the
processing of the data being done in the PC via asp. The routine
will probably not work for the older data taken with the vme array
processors (the record formats/header info was changed).
2. The routine needs to be positioned at the start of a header. If a
file starts with a partial record, use:
rew,lun
searchhdr(lun)
This will position you at the start of the first header. The scan keyword
allows for random positioning within the file. Unfortunatley, it uses
the scannumber keyword in the header and this header element is not
filled in by the processing programs.
3. Be careful how many records you ask for. You need to have enough
memory to hold all of them. Typical record lenghths are:
pwr : 8628 bytes/rec
mracf: 11404 bytes/rec 128 spclen by 21
clp : 154556 byes/rec (64 spclen by 602)
4. I've tested the routine with mracf,clp, and pwr. Record types of:
tpsd, dspc are not yet implemented. Raw data (meteors etc..) using rectype
rawdat will probably work.
-------------------------------------------------------------------------
DATA FORMAT:
A typical record will contain a header followed by data. The data
format will depend on the type of record (pwr,mracf,...) unless the
/raw keyword is used. An example for a power record is:
d.h STRUCT -> HDRPWR Array[1] header
d.tx FLOAT Array[46] tx samples
d.d FLOAT Array[1600] height data
d.calon FLOAT Array[200] cal on
d.caloff FLOAT Array[200] cal off
The same record read using /raw will look like:
d.H STRUCT -> HDRPWR Array[1] header
d.D1 FLOAT Array[2046] all the data
-------------------------------------------------------------------------
HEADERS:
Each record will have a header containing generic and program specific
entries. The generic parts are :
h.std (standard header), h.ri , h.sps.
The program specific headers are:
h.pwr, h.mracf, h.tpsd,h.clp
A description of the generic headers are listed below
IDL> help,d.h,/st
** Structure HDRPWR, 4 tags, length=444, data length=444:
STD STRUCT -> HDRSTD Array[1] standard header
RI STRUCT -> HDRRIV2 Array[1] ri header
PWR STRUCT -> HDRSECPWR Array[1] program header (pwr )
SPS STRUCT -> HDRSECSPSBG Array[1] sps header
-------------------------------------------------------------------------
THE STD HEADER:
** Structure HDRSTD, 26 tags, length=128, data length=128:
HDRMARKER BYTE Array[4] 'hdr_'
HDRLEN LONG 444 header len bytes
RECLEN LONG 8628 reclen bytes
ID BYTE Array[8] prog id: 'pwr','mracf',etc
VERSION BYTE Array[4] version
DATE LONG 2003185 yyyyddd where ddd is daynumber
TIME LONG 70263 time in seconds from midnite ast
EXPNUMBER LONG 0
SCANNUMBER LONG 0
RECNUMBER LONG 0
STSCANTIME LONG 0
SEC STRUCT -> STRSEC Array[1]
GRPNUM LONG 500 first record of scan for this data
GRPTOTRECS LONG 1
GRPCURREC LONG 1
DATATYPE BYTE Array[4]
AZTTD LONG 1156032 az pointing direction in .0001 deg
GRTTD LONG 150000 za greg in .0001 deg units
CHTTD LONG 150002 za ch in .0001 deg units
POSTMMS LONG 70262024 millisec from midnite for positions
-------------------------------------------------------------------------
THE RI HEADER:
** Structure HDRRIV2, 12 tags, length=48, data length=48:
EXTTIMING LONG 1 using sps
SMPMODE LONG 1 use gw pulses
PACKING LONG 12 12 bit sampling
MUXANDSUBCYCLDE LONG 0
FIFONUM LONG 12 use both chan 1,2 (dual beam)
SMPPAIRIPP LONG 2046 total samples 1 ipp
IPPSPERBUF LONG 2 ipps per output rec
IPPNUMSTARTBUF LONG 999 ipp number from start
IPP FLOAT 10000.0 ipp in usec
GW FLOAT 2.00000 gw in usec
STARTON LONG 0
FREE LONG 544368000
note: some of the ri header is duplicated in the sps header.
when in doubt, use the sps header info.
-------------------------------------------------------------------------
THE SPS HEADER:
** Structure HDRSECSPSBG, 16 tags, length=212, data length=212:
ID BYTE Array[4] section id
VER BYTE Array[4] version
IPP FLOAT 10000.0 ipp in usecs
GW FLOAT 2.00000 sample time in usecs
BAUDLEN FLOAT 4.00000 baudlen of code in usecs
BWCODEMHZ FLOAT 0.250000 bandwidth of code in mhz
CODELENUSEC FLOAT 52.0000 codelen in usecs
TXIPPTORFON FLOAT 373.000 tx ipp to rf on in usecs
RFLEN FLOAT 52.0000 rf len in usecs
NUMRFPULSES LONG 1 number of rf pulses
MPUNIT FLOAT 52.0000 multipulse unit
MPSEQ LONG Array[20] multipulse seq (if mpunt gt 1)
CODENAME BYTE Array[20] name of code used
SMPINTXPULSE LONG 46 samples taken in tx pulse
NUMRCVWIN LONG 2 number of receive windows.
RCVWIN STRUCT -> HDRSPSRCVWIN Array[5]
The SPS.RCVWIN structure:
IDL> help,d.h.sps.rcvwin,/st
** Structure HDRSPSRCVWIN, 3 tags, length=12, data length=12:
STARTUSEC FLOAT 300.000 usecs from start of rf on
NUMSAMPLES LONG 1600 number of sample taken
NUMSAMPLESCAL LONG 0 number samples that are cal samples
-------------------------------------------------------------------------
POWER PROFILE RECORDS:
istat=atmget(lun,d,rectype='pwr',nrec=10)
help,d,/st
Structure <829c52c>, 5 tags, length=8628, data length=8628, refs=1:
H STRUCT -> HDRPWR Array[1] header
TX FLOAT Array[46] tx samples
D FLOAT Array[1600] data samples
CALON FLOAT Array[200] cal ON samples
CALOFF FLOAT Array[200] cal off samples
the pwr header d.h.pwr contains:
IDL> help,d.h.pwr,/st
** Structure HDRSECPWR, 14 tags, length=56, data length=56:
ID BYTE Array[4]
VER BYTE Array[4]
PROGMODE LONG 1000
DCDMODE LONG 0
RECTYPE LONG 0
TXSMPSCALE FLOAT 0.00000
SPCRECSPERGRP LONG 0
SPCCURREC LONG 0
HIPPSAVGED LONG 1000
SPCNUMHEIGHT LONG 0
SPC1STHEIGHT LONG 0
SPCLENFFT LONG 0
SPCAVGED LONG 0
SPCTHISREC LONG 0
only the hippsavged (number of ipps averaged) is filled in.
-------------------------------------------------------------------------
MRACF RECORDS:
print,atmget(lun,d,rectype='mracf')
IDL> help,d,/st
* Structure <8266f94>, 5 tags, length=11404, data length=11404, refs=1:
H STRUCT -> HDRMRACF Array[1] header
TXSPC FLOAT Array[128] spectra of tx samples
DSPC FLOAT Array[128, 16] 16 data spectra of 128 pnts
NSPC FLOAT Array[128, 4] 4 noise spectra of 128 pnts
DC FLOAT Array[46] dc points (not returned)
note: the last data spectra d.dspc[*,15] is not computed correctly
(it is overlapped with the first noise spectra).
The mracf header d.h.mracf contains:
IDL> help,d.h.mracf,/st
** Structure HDRSECMRACF, 20 tags, length=80, data length=80:
ID BYTE Array[4]
VER BYTE Array[4]
IPPSAVGDATA LONG 1000 ipps averaged
IPPSAVGNOISE LONG 0
NUMIFFREQ LONG 0
NUMHEIGHTS LONG 20 number of heights
NUMLAGS LONG 64 number of lags. spclen=2*numlags
NUMDCPNTS LONG 0
FIRSTTXSMP LONG 0
RECISDATAREC LONG 0
HEIGHTSTHISREC LONG 0
DCPNTSTHISREC LONG 0
TXSPIPP LONG 0
TXHEIGHT LONG 0
NUMFREQSW LONG 0
TXFRQOFF1 FLOAT 0.00000
TXFRQOFF2 FLOAT 0.00000
TXFRQOFF3 FLOAT 0.00000
TXFRQOFF4 FLOAT 0.00000
FR3 LONG 0
the only elements returned are the 3 listed.
-------------------------------------------------------------------------
CODED LONG PULSE RECORDS:
print,atmget(lun,d,rectype='clp')
IDL> help,d,/st
* Structure <82630e4>, 2 tags, length=154556, data length=154556, refs=1:
H STRUCT -> HDRCLP Array[1] header
DSPC FLOAT Array[64, 602] 602 spectra of len 64
The coded long pulse header contains:
IDL> help,d.h.clp,/st
*Structure HDRSECCLP, 14 tags, length=56, data length=56:
ID BYTE Array[4]
VER BYTE Array[4]
IPPSAVG1FREQ LONG 1000 ipps averaged
NUMIFFREQ LONG 0
NUMHEIGHTS LONG 602 number of heights
SPC1STHEIGHT LONG 0
SPCHEIGHTSTEP LONG 0
DECIMATEFACTOR LONG 8
ZEROEXTDXFORM LONG 0
FIRSTTXSMP LONG 0
SPCLEN LONG 64 length of spectra
DECIMATEDCODELEN
LONG 0
SPCTHISREC LONG 0
FILL LONG 0
note: the current implementation assumes there are no tx or cal
samples returned.
-------------------------------------------------------------------------
RAWDAT RECORDS:
print,atmget(lun,d,rectype='rawdat')
IDL> help,d,/st
* Structure <827d394>, 3 tags, length=123268, data length=123268, refs=1:
H STRUCT -> HDRRD Array[1] header
D1 COMPLEX Array[7680] channel 1 ch if dual beam
D2 COMPLEX Array[7680] channel 2 gr if dual beam
The rawdat header contains:
IDL> help,d.h,/st
* Structure HDRRD, 3 tags, length=388, data length=388:
STD STRUCT -> HDRSTD Array[1]
RI STRUCT -> HDRRIV2 Array[1]
SPS STRUCT -> HDRSECSPSBG Array[1]
The rawdat program always returns the data as a complex array. It does
not split the data up into cal or noise (as if /raw was always set).
(See /pkg/rsi/local/libao/phil/atm/atmget.pro)
ATMHIST - COMPUTE HISTOGRAM OF SAMPLED VOLTAGES
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
atmhist - compute histogram of sampled voltages
SYNTAX: npnts=atmhist(fname,histar,maxtm=maxtm,maxrec=maxrec,maxsmp=maxsmp,$
verb=verb,d=d)
ARGS:
fname: string filename to process (must include complete path)
histar[4096,n]: long contains the returned histograms of the
voltage samples.
n=2 if 1 channel is used,
n=4 if 2 channels are used (ch and dome)
KEYWORDS:
maxtm: float limit histogram to maxtm seconds of data
If nothing specified then maxtm is defaulted to 100 secs.
maxrec: float limit to maxrec records of data
maxsmp: float limit to maxsmp samples of data.
verb : if set then plot the progress in 10% steps..
outputing the ipps and samples read so far.
RETURNS:
npts: double number of counts is a single histogram.
histar[4096,n] long histogram of the voltages.
d[100]: if supplied then return the first 100 records read
in d
DESCRIPTION:
Compute a histogram of the voltages sampled from an aeronomy rawdata
file. Any transmitter samples are discarded. Any cal samples are included.
The number of samples read depends on the keywords maxtm, maxrec, or
maxsmp. If none of these are set then 100 seconds of data are used.
The routine will read 100 records at a time computing a cumulative
histogram for the sampled values. The assumptions made by the
program are:
1. the data is rawdat voltages.
2. All of the records read are the same as the first record of the
file.
3. The ri is set to 12 bit sampling. Actually this is not a real
limitation, just that the returned values will only cover part of the
4096 length histogram.
4. The counts are stored as longs so don't have more than 2^31-1
counts in any bin.
5. The routine alway reads a multiple of 100 recs so you may get a little
more data than you asked for.
The returned data is in histar[4096,n] where:
0:4095 maps into digitizer levels -2048 to +2047
n is:
for 1 channel (dome or ch) n=2
for 2 channels (dome and ch) n=4
For any pair
histar[*,0] is the digitizer sampling the real part of the complex num
q digitizer (right bnc of the pair)
histar[*,1] is the digitizer sampling the img part of the complex num
i digitizer (left bnc of the pair)
EXAMPLES:
file='/share/aeron5/t2163_19mar2006.003'
npts=atmhist(file,histar,maxtm=200,/verb)
x=fingen(4096) - 2048
plot,x,histar[*,0]
oplot,x,histar[*,1],color=colph[2]
(See /pkg/rsi/local/libao/phil/atm/atmhist.pro)
ATMPWRCMP - DECODE/COMPUTE POWER PROFILE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
atmpwrcmp - decode/compute power profile
SYNTAX: istat=atmpwrcmp(d,dcdpwr,usecal=usecal)
ARGS:
d[n] :{rd} array of rawdat power profile data.
KEYWORDS:
usecal: if set then decode and return the cal with the data.
RETURNS:
istat :int 1 ok, 0 error
dcdpwr[m]:float array holding decoded power data starting at first
height (d.h.sps.rcvwin[0].startusec)
DESCRIPTION:
Decode and compute power for power profile data using the barker
code or the 88 baud length code. The decoding is done with the
theoretical codes (not the transmitter samples). All of the
ipps in d[n] are averaged together. By default the routine
decodes only the first receive window. If /usecal is included, then
the cal receive window (if it exists) will be part of the decoding.
EXAMPLE:
usedome=0
istat=atmget(lun,d,nrecs=1000,nrectype='rpwrb')
istat=atmcmppwr(d,dcdpwr,/usecal)
; figure out how to label the plot
za=(keyword_set(usedome))?d[0].h.std.grttd/10000. $
:d[0].h.std.chttd/10000.
lab=string(format=$
'("power profile tm:",a," az:",f6.2," scan:",i9)',$
fisecmidhms3(d[0].h.std.time),d[0].h.std.azttd*.0001,d[0].h.std.scannumber)
; compute the height from the range
range=(findgen(n_elements(dcdpwr))*d[0].h.sps.gw + d[0].h.sps.rcvwin[0].startusec)*.15
hght =range*cos(za*!dtor)
maxhght=(long(hght[n_elements(hght)-1L])+49L)/50L *50L
maxpow=max(dcdpwr)
maxpow=long(maxpow+49)/50l * 50
hor,0,maxpow
ver,0,maxhght
plot,dcdpwr,hght,title=lab,xtitle='power',ytitle='altitude [km]'
(See /pkg/rsi/local/libao/phil/atm/atmpwrcmp.pro)
ATMPWRSCANCAL - GRAB CALS FOR A SET OF RAWDAT POWER SCANS
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
atmpwrscancal - grab cals for a set of rawdat power scans
SYNTAX: n=atmpwrscancal(fbase,fnum1,fnum2,calD,maxblks=maxblks,verb=verb,yscale=yscale)
ARGS:
fbase :string base name for files to scan. includes directory and fname
upto the filenumber: eg:"/share/aeron5/t1193_04sep2012.".
fnum1 : int first filenumber to scan
fnum2 : int last filenumber to scan
KEYWORDS:
maxBlks: long maximum number of blocks to return. The default is 1000.
verb : 1-output filenames and block numbers as they are processed
2-output 1 and plot the cal on,off for each 10 sec block
yver :[min,max] If verb=2 then the the vertical scale for the plotting.
If not provided, the default is [1000,40000]. The plot uses
/ylog so you can see the on,and off.
RETURNS:
n:int > 0 number of entries in calD[]
0 error
-1 - could not open first file fbase+fnum1
-2 - no power data found in first file
calD[n]:{} caldata for each 10sec block
{ fnum: 0L,$
blk : 0l,$ ; in file
h :d.h,$ ; header for first rec of 10 sec block
nipp: 0L,$ ; we found in this block
nfifos:0L,$ ; number of fifos 1,2..
calOn:fltarr(maxipp,nfifos),$ ; holds calOn totPwr each ipp
0=fifo1,1=fifo2. only nipp valid
calOff:fltarr(maxipp,nfifos),$ ; holds calOff totPwr each ipp
0=fifo1,1=fifo2
}
DESCRIPTION:
Scan rawdat files for power profile blocks (10secs) of data. Compute the power in the
calOn,calOff, and return it in the calD[n] structure. Each element of calD[] holds the
information for an entire 10 sec block. A block read will finish after 1000 ipps, or a data record
of a different type (say mracf is found).
For each block of data, compute the total power for the calOn and calOff of each ipp.
Then average over the calOn, and caloff samples (giving 1 calon,1 caloff value for each ipp,
fifoChan). When averaging, ignore 25 samples at the start of calOn,calOff. This excludes the
calOn to calOff transition (which differs for the dome and the ch).
The returned calD[n] array of structures contains:
calD[n].fnum the filenumber where this block came from
calD[n].blk the 10sec block within the file where data came from (count from 0)
calD[n].h the header from the first ipp of the block.
calD[n].nipp the number of valid ipps in calOn,calOff for this block (in case the
block had less than 1000 ipps before a different record type was found.
calD[n].nfifos the number of fifos found. This determines the 2nd dimension of
calOn[x,nfifos],calOff[x,nfifos]
calD[n].calOn[maxipp,nfifos] Hold the calOn total power for each ipp. [x,0] is fifo1 (ch)
[x,1] is fifo2 (dome).
calD[n].calOff[maxipp,nfifos] Hold the calOff total power for each ipp.
Notes:
0. nfifos, and ipps/buf are taken from the first block of the first file. If this
configuration changes (between fnum1,fnum2) then it will still use the values from
the first rec of first block of first file.
1. I've used 10secs of data in the above. It actually is 1000 ipps. If the power ipp is
different, then the time of each block will will be 1000*ipp.
2. If there are 1000 +n ipps,then a different data rec, you will get a second
block with only N ipps (maybe there should be a keyword for specifying the minimum
number of ipps for a block?).
3. The ipps within a block will be contiguous in time (spaced by an ipp). There will
normally be time gaps between blocks (if they run more that just pwr in the cycle).
You can look at the time in the header to see the start of each block.
EXAMPLE:
fbase='/share/aeron5//share/aeron5/t1193_04sep2012.' ; note the . a the end
n=atmpwrscancal(fbase,71,99,calD,/verb)
(See /pkg/rsi/local/libao/phil/atm/atmpwrscancal.pro)
METMON - MONITOR METEOR DATA
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
metmon - monitor meteor data
SYNTAX: img=metmon(fbase,num,lun=lun,toavg=toavg,numDisp=numDisp,$
useCh=useCh,d=d,val=val,pixwin=pixwin,ilim=ilim,$
dirhc=dirhc,hc=hc,invert=invert)
ARGS:
fbase: string Base filename to read from (eg '/share/aeron5/24jul03').
Do not include the . or file number.
num : int file number of file to start on.
KEYWORDS:
lun : int If provided, then process this file from the current
position. When done with the file return. In this case
fbase and num are ignored. You also can not jump to other
files from the internal menu.
toavg: long number of ipps to average together. The default is 5
numDisp: long The number of averaged ipps to display in the image. The
default is 900 averaged ipps.
useCh: If set then start displaying the first data channel. This
is the carriage house in dual beam experiments. The default
for dual beam is the dome. You can change this from the
internal menu.
val[2]: float Clip the power (a/d levels squared) to [min,max] and
then scale to the full range of the display device (
0 to 255). The default is to clip the power at 6 sigma
(as measured from the first image).
pixwin: if set then use a pixwin when drawing a new image. This
cuts down on the flashing. It is useful when you are
averaging only 1 or 2 ipps.
ilim : long if supplied then limit image to these indices of the
data rec (0 based).
dirhc : string directory to start hardcopy
hc : in 1 hardcopy on, 0 hardcopy off
invert: if set then reverse black, white in lut.
RETURNS:
img[m,numDisp]:float the last image displayed.
d[n] if this keyword is provided then the data read
for the last image is returned here.
DESCRIPTION:
metmon will make continuous images of meteor data. The mean is removed from
the voltages, power is computed, and then toavg ipps are averaged together.
numDisp averaged ipps are then displayed as height versus time. When
the end of file is hit, the routine will advance to the next filenumber.
If no data is available, the routine will wait until it becomes available.
The user can modify things by hitting any key and bringing up an
internal menu:
Cmd CurVal function
a dome antenna ch, dome
f 70 move to new fileNum
h 0,1 hardcopy on,off
hd dirName directory name for hardcopy
l list all files
n next file (or quit if 1 file)
p recnum position to rec in file
pr 0,1 display profiles in image
q to quit
r rewind current file
s 0 single step 0,1 (off,on)
sc 6 rescale images to nsigmas
cur trac cursor
d debug. stop in metmon. .continue to continue
otherkey continue
The commands are:
a dome/ch This lets you switch between dome and carriage house display
f filenum You can start displaying at the start of a different
filenumber. If the new filenumber is illegal then no
change is made.
l This will list all of the available files starting with
fbase. The last file will also contain its size.
n move to the next file number.
p recnum position to rec in file. Rec num counts from 1
pr 0,1 Use cursor to display profiles in image. 1 on, 0 off.
q quit the routine.
r rewind and start over in the current file.
s 0,1 turn on,off single step mode.When it is on, the routine
will pause after every image waiting for the user to
enter return.
d debug. This will stop in the metmon routine.
otherKey any other key will cause the display to continue.
This allows you to pause the display to look at it for
a while.
EXAMPLES:
An example of using this routine at AO:
1. slogin to either fusion00,01, or 02
2. idl
3. @phil
4. @atminit
5. xloadct .. then click on bw linear for scale color table
6. setup the parameters for the call..
file='/share/aeron5/23Dec03T1748' .. this should be the path and
base filename for the files to use.
usech=1 .. set this to 0 if you want the dome output
num=272 .. first file number of your experiment to use.
7.start the program:
img=metmon(file,num,usech=usech)
8. Hit any key to bring up the keyboard menu.
Possible problems:
By default the lut (color lookup table) is scaled to 6 sigma of the
first image displayed. If this image has a large meteor in it, the
other plots may not have the correct contrast. You can fix this by:
a. hit any key;
b sc nsig .. where nsig is the numbers of sigmas to scale the image.
.. it will the compute the new clipping values from the
current image you are looking at.
(See /pkg/rsi/local/libao/phil/atm/metmon.pro)
PLASMACUTOFF - PROCESS SINGLE PULSE PLASMA CUTTOFF DATA
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
plasmacutoff - process single pulse plasma cuttoff data
SYNTAX: istat=plasmacuttof(lun,spc,spcH,spcN,spcLen=spcLen,toavg=toavg,$
freq=freq,flip=flip,nonoise=nonoise,hdr=hdr)
ARGS:
lun: int file descriptor of file to read
KEYWORDS:
spcLen: long length of transform to do. Default is 2048.
there must be at least this many samples of data and/or noise
toavg: long number of ipps to average default is 10000
flip: if set then flip the frequency order of the spectra on return.
nonoise: if set then do not process the noise samples (even if they are
present).
RETURNS:
istat : int 1 returned ok
0 hit eof
lt 0 trouble reading a record (bad header, etc)
spc[m,n,l]: float height spectra with noise subtraction.
m=length of spectra
n=number of height spectra computed
l=number of antennas: 1 if 1 antenna,
2 if both antennas (dome is 2nd)
spcH[m,n,l]: float height spectra without noise subtraction.
spcN[m,l] : float noise spectra
freq[m] : float frequency of the spectra in Mhz.
hdr : {hdr} header from first record read
DESCRIPTION:
Input and process single pulse plasma cutoff data. Each record should
contain only 1 ipp of data. The windows in the ipp can be:
txSamples - skipped
hghtSamples - must be at least spcLen samples. The routine will compute
hghtSamples/spcLen height Spectra
noiseSamples- if present (and nonoise is not set) then there must be
at least spclen noise samples. The routine will compute
(and average together) noiseSamples/spclen noise spectra
for each ipp.
The routine will read toavg ipps. For each ipp it will compute the
height spectra and an average noise spectra. It then averages these
spectra over toavg ipp's.
The height spectra with noise removal is returned in:
spc[spcLen,numHghtSpc,numAntennas].
The height spectra with no noise removed is returned in:
spcH[spcLen,numHghtSpc,numAntennas].
The noise spectra are returned in :
spcN[spcLen,numAntennas].
If the sps buffer has two receive windows, then the 2nd window is taken
as the noise samples. If only one receive window is present, no noise
subtraction is done.
The flip keyword will flip the frequency order of all of the spectra.
If BBM sine (top) goes to digitizer (left) then you probably need
to set /flip.
EXAMPLE:
idl
@phil
@atminit
openr,lun,'/share/aeron5/24Jul03.070',/get_lun
.. avg 10000 ipps
istat=plasmacutoff(lun,spc,spch,spcN,spclen=2048,freq=freq,/flip)
ver,0,5e8
hor
; plot the ch spectra
stripsxy,freq,spcar[*,*,0],0,1e8,/step
; over plot the dome spectra
stripsxy,freq,spcar[*,*,1],0,1e8,/step,/over
NOTE:
The routine will not work:
1. if there are more than 1 ipp per record.
2. if there are fewer than spclen data or noise samples
(See /pkg/rsi/local/libao/phil/atm/plasmacutoff.pro)
SHSCLOSE - CLOSE AN ATM SHS FILE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclose - close an atm shs file
SYNTAX: shsclosel,desc,all=all
ARGS:
desc: {} descriptor returned by shsopen()
KEYWORDS:
all: if set then close all open descriptors.
DESCRIPTION:
close an atm shs file. This frees up the lun.
(See /pkg/rsi/local/libao/phil/atm/shsclose.pro)
SHSCLP - CODED LONG PULSE (RI)
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclp - coded long pulse (ri)
SYNTAX: istat=shsclp(desc,spcBuf1,spcBuf2,ippToAvg=ippToAvg,baudLen=baudLen,
spclen=spclen,txSmpSkip=txSmpSkip,nheights=nheights,$
dinfo=dinfo,hdr=hdr,dotm=dotm,tmi=tmi,verbose=verbose,$
clipRfi=clipRfi,,minipp=minipp,usenoise=usenoise,$
dinp=dinp,firstheight=firstheight,doband=doband)
ARGS:
desc: struct desc that points at data file
KEYWORDS:
ippToAvg: long number of ipps toavg. there are 100 ipps/rec so
should be a multiple of 100 ipps
minipp : long minimum number contiguous ipps to avg, if not found
continue through the file for the next block.
baudlen: float This is actaully the step in the height processing.
It defaults to 1 usec
lets you override it. The units are usecs.
spclen : long The length of the spectra to do. By default it is
rounded up to the next power of 2.
txSmpSkip:float The number of usecs to skip before taking the first
tx sample. This takes in account the filter delay for the
tx samples. The default is 5 usecs.
nheights: long limit to 1..nheights.. default is all
dotm : if set then do detailed tming.. if not, just do
total times.
verbose : if set then output start of each rec and start of
each ipp modulo 10
cliprfi : float if supplied, then clip the spectra of each rfi to
nsig=cliprfi, 0 these voltage points and back xform
(unless usenoise is set)
usenoise : int if true then replace large rms with noise rather
than 0 it.
dinp[m]: int use this at data to process. ippToAvg will be
set to the min of ipptoAvg m*dinp.nipps
firsthght: long if provided then first height to process. count from 0.
default is 0
navgipp : long number of navgipp to return. default is 1.
use this to return unaveraged ipp
set ipptoavg=1 set navgipp=100.. this will return 1 rec
doband : int 1,2, or 12 to do both bands. does not merge bands def:12
RETURNS:
spcbuf1[spclen,nhghts]: float the averaged spectra vs heights for the
first channel (normally ch1 for dual beam).
spcbuf2[spclen,nhghts]: float the averaged spectra vs height if two
channels were used (normally gr for dual beam).
dinfo : {} Structure holding info that was used in the computation
(see below)
dhdr : {} header from the first record averaged.
tmI : {} Timing info (see timing info below for a description). If
dotm=1 then you get detailed info. If not, then just the totals.
DESCRIPTION:
Input and process coded long pulse data taken in raw data mode with
ryans machine. It will read ippToAvg/100 records (default 10) starting
from the current location on disc (pointed to by desc). It
requested number of records into a buffer (so don't make it too large). It currently
does not check if this is a clp record or not (to be added).
For each ipp the processing step is:
1. extract the tx samples using txSmpSkip to determine which
tx Sample to start on.
2. compute the complex conjugate of these samples.
3. for each height:
- multiply the data by the code.
- zero extend to spclen
- fft the data
- compute power and accumulate
- compute how many samples to skip to get to the next starting
height (usually baudlen/sampLen).
4. After accumulating divide by the number of accumulations;
5. Shift each spectra so that dc is in the center (for a 4096 length
xform, shift right by 2048)
6. fill in the dinfo structure with the info that was used for the
computation. Dinfo contaings:
dinfo={ $
gwUsec : gwUsec,$
baudLenUsec: baudLenUsec,$
codeLenGw : codeLenGw,$
nhghts : nhghts,$
hghtStepGw : hghtStep,$
ippAvged : ippAvgd,$
spcLen : spclen,$
txSmpSkip : txSmpSkipL,$
numChn : numChn ,$
bwMhz : bw,$
nsigclip : nsigclip,$
ippI : ippI[ipptoAVg] } ; info on each ipp
ippI hold info on each ipp:
clipI={$
mean:0.,$; spectral mean
rms :0.,$; spectum rms
ndel:0. $, ;nchan removed
rmsLoop:0$; 1 or 2 .. how many times we looped rms
}
a={ rec: 0L,$ in file count from 0
ippRec:0L,$; ipp within rec count from 0
ippCum:0L, ; from start of file count from 0
; following computes on spc.. if nsigdel 0.
clipI : replicate(clipI,2) $; the 2 chan
}
7. also return the header from the first ipp used in the variable hdr.
Interesting info is in the sps portion of the header (although the
baudlen may be incorrect).
TIMING INFO:
Timing info can be returned if tmi=tmi keyword is provided. The tmI
structure holds times for different parts of the code. Each time
structure has the number of times the code was timed, the total sum
and the total sum of squares. The routine printtmI,tmi can be used
to print this info out.
An example output is:
IDL> printtmi,tmi
tmTot:675.096 read: 3.987
CODECONJ Ntot: 1000 tmTot: 0.4328 avgTm0.000433 sig:0.000006
CODEM1 Ntot:601000 tmTot: 51.1495 avgTm0.000085 sig:0.000005
FFT1 Ntot:601000 tmTot:122.8049 avgTm0.000204 sig:0.000008
ACCUM1 Ntot:601000 tmTot:135.5683 avgTm0.000226 sig:0.000006
CODEM2 Ntot:601000 tmTot: 51.2731 avgTm0.000085 sig:0.000004
FFT2 Ntot:601000 tmTot:122.8776 avgTm0.000204 sig:0.000007
ACCUM2 Ntot:601000 tmTot:135.2315 avgTm0.000225 sig:0.000006
BUF1TOT Ntot:601000 tmTot:322.3068 avgTm0.000536 sig:0.000013
BUF2TOT Ntot:601000 tmTot:322.0795 avgTm0.000536 sig:0.000012
All times are in seconds. the columns are:
totmeasTm total wall time
read: total time for reading data.
sectionNm The section of code that was timed.
Ntot is how many times this section was timed.,
tmTot is the total time for this section.
avgTm is the average time for 1 execution of this code
sig is the sigma for the Ntot measurements
The sections are:
CODECONJ : extract and conjugate the tx samples
CODEM1 : extract hght data, multiply by code, 0 extend.
FFT1 : compute the fft (using fttw)
accum1 : compute the power and accumulate
BUF1TOT : total time for buf1. includes some data manipulation times
not included in the above times
The xxx2 : same as xxx1 but for the second chan (if present)
totmeasTm : this is just the sum of the measured times. It does not
include the overhead of timing (about 2usecs per call).
The example above had 1000 codeconj so 1000 ipps were processed (10 secs
of data).
You'll notice that computing the power and accumulating is taking
longer than the fft.
For this computer fftw in C took about 44 usecs for a 1k xform. So it
should take about 211 usecs for a 4k. It is averaging 204 so there doesn't
look like there is any speed up time in the fft.
Code speedup would probably come by writing an external module
that was passed a tx rec, hght rec, and then it did the fftw and
power, accumulate, returing the averaged buf. This might give you
50% speed up.
(See /pkg/rsi/local/libao/phil/atm/shsclp.pro)
SHSCLPCHKIPP - SEE IF THIS IS CLP IPP
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclpchkipp - see if this is clp ipp
SYNTAX: istat=shsclpchkipp(txbuf,smpUsec=smpUsec,codeUsec=codeUsec,$
corRatio=corRatio,txminpwr=txminpwr)
ARGS:
txbuf[2,N]: short tx buffer i,q n samples
KEYWORDS:
smpUsec : float sample rate in usecs (def=.2 usecs)
codeUsec : float codelength in usecs (def=440 usecs)
RETURNS:
istat[n] : 1 - is a ipp
0 not a clp ipp
one for each tx buf passed in
corRatio : float acf of txbuf . 0..1
DESCRIPTION:
Check to see if this ipp is a clp ipp by looking at the power across the
transmitter pulse
values used:
midipp: 150.. must have power .. gets rid of power profile which may have
echo at 440,500
power : 52 usecs
clpipp: 440.
topSid: 500 usecs
mracf : 300 usecs
txSmpDelayUsecs .. usually about 5.6 usecs
slopEnd : 5 usecs; .. move this far back from falling edge of pulse before average start
avgSmpls:100= 20 secs
make sure power at 440 - 6 usecs
acf:
- compute acf for lags 0 .. 50 (0 to 10 usecs)
- avg lags 20-50 ( 4 to 10 usecs)
- normalize to 0 lag
make sure correlation is less than acfLim (.4)
- this gets rid of topsid.. so i'll remove the topside limits
- mracf, pwr, clp all have low correlation
(See /pkg/rsi/local/libao/phil/atm/shsclpchkipp.pro)
SHSCLPIMGMONSCAN - SCAN FOR PROCESSED CLP FILES
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclpimgmonscan - scan for processed clp files
SYNTAX: nfiles=shsclpimgmonscan(dir,fbase,fiar,firsttime=firsttime,
allhdrs=allhdrs,maximgfile=maximgfile
ARGS:
dir : string to search for .hdr files. Include the
'/' at the end of the dir name
fbase : string for filename to look for. Include up to the
_ before the filenumber digits (see below)
KEYWORDS:
nrec: long number of records to process (default is 1000).
warning.. this is the number of records, not ipps
unless ipps/buf = 1.
firstTime: if set then this is the first time the routine
is called with this dir,fbase. If not set then
fiar will have the info from the previous call.
the routine will add new info to this array.
allhdrs: if set then return all header for each file.
Each .shs file will have multiple 10 sec
images generated from it. By default only the
first hdr of the .shs file is returned in fiar[i]
With allhdrs set, all the headers are returned
in fiar[i].hdrs[20]. Leaving this off speeds up
the scanning. If you need the az,za,time exactly
for each img then turn this on. It only works
with firsttime=1
maximgfile:long max number of images we can have in 1 input file
def: 120 (or 1 sec integrations)
RETURNS:
nfiles: long number of raw input files found
fiar[nfiles]:{} stuct holding info on all the image files found:
DESCRIPTION:
Scan for .shs processed image files. The routine looks for
the .hdr files in the provided dir, fbase. It returns the number
of input files found (note that there are multipled 10sec img files
for each raw input file.
EXAMPLES:
dir='/raid/radar/phil/t1193_20120120/'
fbase='t1193_20jan2012_'
n=shsclpimgmonscan(dir,fbase,fiar,first=1)
Structure <170f688>, 4 tags, length=196, data length=196, refs=1:
FNUM LONG 0
NIMG LONG 10
HDR STRUCT -> Array[1]
IMGNUMAR LONG Array[20]
n=shsclpimgmonscan(dir,fbase,fiar,first=1,/allhdrs)
help,fiar,/st
Structure <19c57b8>, 4 tags, length=2248, data length=2248, refs=1:
p FNUM LONG 0
NIMG LONG 10
HDR STRUCT -> Array[20]
IMGNUMAR LONG Array[20]
Note that without /allhdrs, you only get the header from the
first 10 sec block of the file.
(See /pkg/rsi/local/libao/phil/atm/shsclpimgmonscan.pro)
SHSCLPPARSETMING - PARSE TIMING INFO FROM CLP HDR FILES
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclpparsetming - parse timing info from clp hdr files
SYNTAX:n=shsclpparsetming(dirName,tmIAr)
ARGS:
dirName: string name of directory holding .hdr files
RETURNS:
n: long number of entries found
tmIAr[n]: {} holds timing info:
DESCRIPTION:
The clp processing generates a .dcd and a .hdr file for each
integration (usually 10 seconds). This routine will search
all of the .hdr files in a directory and return the timing
info. This info tells how long it took to complete different
parts of the processing.
The tmIAr struct has:
;tmI={ wstS: 0. , $; secs waited before start
fftU: 0. , $; usecs for 1 fft
accumU: 0. , $; usecs for power and accum 1 height,1 ipp
ippS: 0. , $; secs for 1 ipp
inpS: 0. , $; secs for input
outS: 0. , $; write out time secs
tmTot:0.} ; tot tm secs this block
(See /pkg/rsi/local/libao/phil/atm/shsclpparsetming.pro)
SHSCLPQUICKLOOK - COMPUTE A CLP SPECTRA AND DISPLAY IMAGE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsclpquicklook - compute a clp spectra and display image
SYNTAX: shsclpquicklook,fnum,fbase,dir0=dir0,$
fband=fband,doband=doband,
firstHghtKm=firsthghtKm,lastHghtKm=lastHghtKm,hghtStepUsec=hghtStepUsec,$
zx=zx,zy=zy,$
dinfo=dinfo,hdr=hdr,img1=img1,img2=img2 ,
ARGS:
fnum: long file number to process. -1 --> last file
fbase: string name of subdir in dir0 to search. eg let fbase='expabc'
filenum =-1 searches dir0+'expabc/*.shs
filenum =val searches dir0+'expabc/expabc_nnn.shs'
KEYWORDS:
dir0: string base dir. default: /net/rserv/export/RSERVLV2/daeron/'
fband: float[2] frequenciy offsets (in Mhz) for bands. def:[-5.5,5.5]
doband: int bands to process. 1->band 1, 2-> band 2, 12-> band 1 and 2
(12 does not currently merge images).
firstHghtKm:float first height to display. def=90km.
If less than first height then first height will be used.
lastHghtKm:float last height to display. def=450Km
If greater than last height recorded, then clip to last height.
hghtStepUsec:float Spacing between heights (in usec, not km). def: 2usec
the baud is normally 2usec.
zx :float zoom factor in x. def: -5 if zx=0 then use default
zy :float zoom factor in y. def: search for zy giving le 790 pixels
RETURNS:
dinfo: {} struct holding info used in compuatation
hdr : {} header read from file
img1[n,m]: float img 1 (if selected). this has not been zoomed.
img2[n,m]: float img 2 (if selected). this has not been zoomed.
DESCRIPTION:
Read and process an echotek clp image. It tries to read
10 secs of data. It will always start at the begining of the
selected file.
To select the file:
- dir0: is the directory holding the experiment directories
the default is /net/rserv/export/RSERVLV2/daeron/. This the destination
for the copy from daeron to aonet.
Another alternative would be to read it directly from daeron:
/net/daeron/data/
- fbase: is the base filename
this is the subdirectory in dir0 and well as the start of the filename.
eg: 'hf_20180315'
- fnum: this is the file number in the subdirectory to input
fnum=-1 will use the most recent file.
Bands to process:
There are normally 2 bands [-5.5,5.5] if up,down or [5.5,9.5] upper band
The parameter doband lets you select which bands to process.
doband=1 --> process only band 1, doband=2 --> process only band 2,
doband=12 --> process both bands. By default both bands are processed.
Processing fewer bands speeds up the processing.
Heights to process
By default all heights 90 through 450 km are processed with 2 usec spacing (same as the def baud)
firstHghtKm specifies a different first height, Asking for a height before the first available
height was def to the first height.
lastHghtKm specifies a different ending height. asking for a height greater than available will
default to the last available height
heightStepUsec: spacing between heights.
image zoom
zx,zy - x,y image zoom. Positive make the image bigger, neg shrinks the dimension (zoom needs to be integral).
In x there are 4096 channels. a zoom of -4 would make an image with 1024 horizontal pixels.
In y, a standard can support 795 pixels (the axis are extra).
if zx=0 or is not supplied then it defaults to -4.
if zy=0 or is not supplied then the zy is computed such that there are as many pixels as possible and does
not exceed 795 pixels.
eg with 90km to 450 km there are (450-90)/.15(km/usec) = 2400 pixels. It would compute zy=-4 giving 600 pixels.
Returned info.
dinfo
IDL> help,dinfo,/st
** Structure <1c29b08>, 13 tags, length=44056, data length=40050, refs=1:
GWUSEC DOUBLE 0.20000000
BAUDLENUSEC FLOAT 2.00000
CODELENGW LONG 2200
NHGHTS LONG 1201
HGHTSTEPGW LONG 10
IPPAVGED LONG 1000
SPCLEN INT 4096
TXSMPSKIP LONG 25
NUMCHN INT 2
BWMHZ DOUBLE 5.0000000
NSIGCLIP FLOAT 0.00000
USENOISE INT 0
IPPI STRUCT -> Array[1000]
hdr: {}
IDL> help,hdr,/st
** Structure <21b62f8>, 18 tags, length=104, data length=90, refs=1:
TABLEREC LONG 0
TABLESIZE LONG 9844000
DATAWIDTH INT 2
DATATYPE STRING 'short'
NUMDIMS INT 2
DIM0 LONG 49220
DIM1 LONG 200
NUMCHANNELS INT 2
TXSTART LONG 0
TXLEN LONG 5200
DATASTART LONG 6000
DATALEN LONG 44000
NOISESTART LONG 50000
NOISELEN LONG 20
SYSTIME STRING '09:30:07'
AZ FLOAT 347.997
ZAGREG FLOAT 1.10000
ZACH FLOAT 12.0000
img1[n,m] : float image band 1 (if processed)
img2[n,m] : float image band 2 (if processed)
Time to compute:
cpu time nbands hgtStep hghtrange
min usec km
megs3 3:30 2 4 90 -690
megs3 2:10 2 4 90-450
megs3 1:07 1 4 90-450
megs3 2:12 1 2 90-450
(See /pkg/rsi/local/libao/phil/atm/shsclpquicklook.pro)
SHSCMPREC - COMPUTE CURRENT RECORD POSITION
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shscmprec - compute current record position
SYNTAX: currec=shscmprec(desc)
RETURNS:
currentrec: long current record we are about to read.
count from 0
DESCRIPTION:
Compute the record we are about to read. count from
0. If we are in the middle of a record, return the current record
(removing the fractional part
from the first record of the file.
(See /pkg/rsi/local/libao/phil/atm/shscmprec.pro)
SHSFINDTM - FIND FILE,REC FOR A GIVEN TIME
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsfindtm - find file,rec for a given time
SYNTAX: istat=shsfindtm(yymmdd,hhmmss,fI,indFile,indrec)
ARGS:
yymmdd: long date we want
hhmmss: long time (ast) we want
fI[n] : {} array of struct returned from shsscandir()
RETURNS:
istat : int 0 this date does not occur in dir
-1 time is before beginning of first file
-2 time is after the end of last file
indFile: long index into fI[] for file we want.. cnts from 0
INDREC : long record we want in file (cnt from 0)
(See /pkg/rsi/local/libao/phil/atm/shsfindtm.pro)
SHSGET - INPUT ATM SHS DATA
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsget - input atm shs data
SYNTAX: istat=shsget(desc,d,nrec=nrec,posrec=posrec)
ARGS:
desc: {} stucture returned by shsopen
KEYWORDS:
nrec: long number of records (hdr/data) to input
posrec: long position to this record before reading (count from 0).
if posrec is not supplied or it is set < 0 then the
next record will be read.
RETURNS:
istat: > 0 number of recs found
d[nrec]: {shsdat} structure holding data
DESCRIPTION:
Input data from the file pointed to by desc (desc is a structure
returned by shsopen()). The data is returned in the structure d. d will
hold one record from the .shs file. A record will have multiple ipps
(for 10 millisec ipps there are 100 records .. 1 second of data).
If nrec=n is provided with n gt 1 then d[n] will be an array.
The user can optionally position in the file before reading
with the posrec= keywword.
The data format on disc is:
primary header
datahdr
data
slop
EXAMPLES:
file='/mnt/daeron/testCont/meteor_08oct2011_000.shs'
istat=shsopen(file,desc)
; read in the next record:
nrec=shsget(lun,d)
; The d struct contains:
help,d,/st
DHDR STRUCT -> Array[1]
NIPPS LONG 100
D1 STRUCT -> Array[1]
D2 STRUCT -> Array[1]
help,d.d1,/st
TX INT Array[2, 100, 100]
DAT INT Array[2, 16100, 100]
CAL INT Array[2, 10, 100]
D1 holds the data for channel 1, D2 has the data for channel 2.
The data is returned as complex so:
d1.dat[i0,i1,i2] where
i0 - 2= the i,q complex samples
i1 - 16100= the data samples in 1 ipp
i3 - 100 = the 100 ipps in a record
To read multiple records:
nrec=shsget(lun,d,nrec=3)
; in this case d will be an array of structs: d[3]...
NOTES:
31aug11.. i was interpreting the tx skip variables incorrectly
the txskip, only tell you
the timing offsets. the data is packed together
txlen,datalen,noiseLen
(See /pkg/rsi/local/libao/phil/atm/shsget.pro)
SHSGETSPC - INPUT DECODED SPECTRA
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsgetspc - input decoded spectra
SYNTAX: nspc=shsgetspc(dir,fbase,spcar,freq,hdr)
ARGS:
dir :string directory holding file
fbase :string file basename (see below)
RETURNS:
nspc : long < 0 error..
> 0 number of spectra in array
spcar[nchn,nspc]:float spectra
freq[nchn]: float freq for spectra
hdr : {} header info
DESCRIPTION:
Read a shs decoded averaged spectra file (typically avged for 10 secs).
This has been processed by clp1shs.c program.
The user inputs the directory and base filename (see below).
The spectra, freq array, and header struct will be returned.
This routine is normally used for debugging...
EXAMPLE:
Suppose the integrated files you want are:
/net/rserv2/data/atm/clp/test_20170101_001.hdr
/net/rserv2/data/atm/clp/test_20170101_001.dcd
The call would be:
dir ='/net/rserv2/data/atm/clp/'
fbase='test_20170101_001' ; note no trailing .
nspc=shsgetspc(dir,fbase,spcar,freq,hdr)
(See /pkg/rsi/local/libao/phil/atm/shsgetspc.pro)
SHSHDR_DAT - GET DATA HEADER FOR TABLE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shshdr_dat - get data header for table
SYNTAX: istat=shshdr_dat(lun,dhdr,ahdr=ahdr)
ARGS:
lun: int lun assigned to file vis shsopen. desc.lun
RETURNS:
istat: int 1 - got header
0 - i/o error or eof
-1 - never saw start of header "HEADER:"
-2 - did not find end of header "]"
dhdr: {} strucuture holding primary header
adhr[]: string array holding ascii lines input
DESCRIPTION:
Read in the primary header from the file pointed to by desc.lun. The
file will be rewound before reading is started. The file will be
left positioned after the start of the primary header.
(See /pkg/rsi/local/libao/phil/atm/shshdr_dat.pro)
SHSHDR_PRI - GET PRIMARY HEADER OF FILE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shshdr_pri - get primary header of file
SYNTAX: istat=shshdr_pri(lun,phdr,ahdr=ahdr)
ARGS:
lun: int file descriptor pointing to data file
RETURNS:
phdr: {} strucuture holding primary header
adhr[]: string array holding ascii lines input
DESCRIPTION:
Read in the primary header from the file pointed to by lun. The
file will be rewound before reading is started. The file will be
left positioned after the start of the primary header.
(See /pkg/rsi/local/libao/phil/atm/shshdr_pri.pro)
SHSMERGEIMG - MERGE TWO IMAGES AVERGING OVERLAP
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsmergeimg - merge two images averging overlap
SYNTAX: istat=shsmergeimg(img1,img2,cfrAr,bw,imgMerged,imgI,nedgslop=nedgeslop)
ARGS:
img1[n,m]: float one img to merge
img2[n,m]: float 2nd img to merge
cfrar[2]: double freq, Mhz center each band.
if 4096 chan, the freq for chan[2048] cnt from 0
bw : double bw in mhz for each image
KEYWORDS:
edgslop: long number of channels to ignore from each edge of each
img. def is 34. this area has no data.
RETURNS:
istat : 0 ok
-1 error
imgMerged[nn,m]: float merged image
imgI : {} holds info about the new img
- imgI.nchan
- imgI.f0 first freq
- imgI.f1 last freq chan
- imgI.df channel freq spacing
DESCRIPION:
merge two images together. It will discard edgeslop channels
from the edges of each input image (def:34).
It will then averge the overlap region and return the final image.
The routine does not reinterpolate the data. If the frequency
step does not divide into the frequency separation of the two bands
then there will be a freq error in the overlap region. this will be a
fraction of the channel width.
The freq returned are:
imgI.f0 - is the first frequency kept from the first image
imgI.f1 - is the last freq kept from the 2nd image
imgI.df - is the original freq freq step
(See /pkg/rsi/local/libao/phil/atm/shsmergeimg.pro)
SHSOPEN - OPEN AN ATM SHS FILE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shsopen - open an atm shs file
SYNTAX: istat=shsopen(filename,desc)
ARGS:
filename: string filename of file to open
RETURNS:
istat: int 1 open ok
0 couldn't open file
desc: {} descriptor holding file info. pass this to the shs io routines.
DESCRIPTION:
Open an atm shs file. Return a file descriptor that you will pass
to the shsxxx i/o routines. The file is left positioned at the start hdr
of the first table.
(See /pkg/rsi/local/libao/phil/atm/shsopen.pro)
SHSPOS - POSITION TO A RECORD IN THE SHS FILE
[Previous Routine]
[Next Routine]
[List of Routines]
NAME:
shspos - position to a record in the shs file
SYNTAX: istat=shspos(desc,rec,currec=currec,nextrec=nextrec)
ARGS:
desc:{} structure returned by shsopen().
rec: long record to position to (count from 0)
KEYWORDS:
currec: if set then reposition to start of current rec
nextrec: if set then position to next record.
RETURNS:
istat: 1 - positioned ok.
0 - beyond end of file
-1 - error positioning
DESCRIPTION:
Position to a record in an shs file. This uses the record length
from the first record of the file.
(See /pkg/rsi/local/libao/phil/atm/shspos.pro)
SHSSCANDIR - SCAN ALL FILES IN AN SHS DIRECTORY
[Previous Routine]
[List of Routines]
NAME:
shsscandir - scan all files in an shs directory
SYNTAX: n=shsscandir(path,fI)
ARGS:
path: string directory path, filename (with wildcards) to search
this string is passed to file_search()
RETURNS:
n : long number of files found
fI[n]: {} holds info for each file
The fI struct contains:
IDL> help,fi,/st
** Structure <8c5cd8>, 5 tags, length=48, data length=44, refs=1:
YYYYMMDD LONG 20150619
FNAME STRING '/net/daeron/data/t1193_20150619b/t1193'...
LTIME STRING '23:51:24'
SECMIDST LONG 85884
NRECS LONG 109
(See /pkg/rsi/local/libao/phil/atm/shsscandir.pro)