NAME:
0atmIntro - Intro to atm routines.
There are two sets of idl routines to process atm data:
- atmxxx To process radar interface taken in the rawdata taking mode.
Use atminit() to initialize idl.
- shsxxx To input/process the echotek digital receiver data.
Use shsinit() to initialize idl.
Starting idl.
- To use the atm routines you need to tell idl where to
look to get these procedures. You can do it manually each
time you run idl, or you can put it in an idl startup file.
Manually:
idl
@phil
@atminit (to process the ri data)
or
@shsinit (to process the echotek data)
Using a idl setup file:
Suppose your home directory is ~jones.
create the file ~jones/.idlstartup
add the line
@phil
to this file.
In your .cshrc file (if you run csh) add the line
setenv IDL_STARTUP ~/.idlstartup
You can then run idl with :
idl
@atminit or @shsinit
You can also put any other commands in the startup
file that you want to be executed each time you
start idl. My startup file contains:
!EDIT_INPUT=500
@phil
Note: if you are running the idl routines at a different site,
replace @phil with the directory at your site.
0. The atmxxx radar interface routines include:
- @atminit : call once after starting idl to initialize the atmxxx
routines.
- openr : use the idl routine open to first open the file.
- atmget : input one or more data records
- atmpwrcmp: decode and compute power for a barker coded power profile
record (you first input the data with atmget().
- atmpwrcmp: decode and compute power for a barker coded power profile
- plasmacutoff: process single pulse plasma line cutoff records.
This routine inputs and processes multiple records.
- metmon : online monitoring of meteor datataking. It can also be
used for offline files.
1. The shsxxx echotek digitial receiver data:
- @shsinit : call once after starting idl to initialize the shhxxx
routines.
- shsopen : open a file to process.
- shspos : position to a record in a file
- shsget : read in 1 or more records
- shsclose : close the file when done.
- Support routines not normally called directly:
- shshdr_pri: input the primary header
- shshdr_dat: input the data header
EXAMPLES using the shs routines:
The file info is kept in a descriptor structure (below it is called
desc). It is generated by the open call shsopen() and should be
freed up with the shsclose call.
A data file contains a primary data header followed by one or
more records. The primary header and data headers are in ascii. The
data is in binary (typically signed shorts).
Records are "tables" in the shs header jargon. A record contains a
data header followed by N ipps of data (typically 100).
To open, read, position, read, close the file you would:
--> idl .. start idl
--> @phil .. if you don't have it in your idl startupfile
--> @shsinit .. init the shs routines (just needs to be called once).
--> file='/share/xserve10/nearField/nearField1/nearField_000.shs'
--> istat=shsopen(file,desc)
; Contents of the descriptor:
help,desc,/st
LUN LONG 100 ; file descriptor from idl openr
FNAME STRING '/share/xserve10/nearField/nearField1/nearF'... filename
FILESIZE LONG 2136692736 ; total number of bytes in file
NUMREC LONG 117 ; total number of records in file
PHDR STRUCT Array[1] ; structure containing primary header
DHDR1 STRUCT Array[1] ; structure containing data header from 1st rec
TBLST LONG 3168 ; byteoffset for start of 1st record
TBLLEN LONG 18262304 ; Length of 1 record (header and data).
; primary header:
help,phdr,/st
VERSION STRING '0.02'
LENDIAN INT 1
DATE STRING 'Aug 11,2005'
TIME STRING '03:20:50 PM AST'
OBSERVER STRING 'MIT'
OBJECT STRING 'PLASMA TEST'
TELESCOPE STRING '430 MHZ'
BWMODE STRING 'ULTRAWIDEBAND'
CHUSED STRING '1 CHANNEL'
COLMODE STRING 'BURST MODE'
GATESRC STRING 'EXTERNAL'
GATEPLR STRING 'POSITIVE'
BURSTCNT LONG 50176
DELAYCNT LONG 0
SYNC STRING 'COHERENT-SIA'
DATAFORMAT STRING '16 BIT COMPLEX'
CLOCKRATE DOUBLE 80.000000
DECIMATE INT 8
OUTPUTRATE DOUBLE 10.000000
RESAMPLER INT 0
IPP LONG 10
SAMPLETIME DOUBLE 0.10000000
; data header (beginning of each record)
help,desc.dhdr1,/st
TABLEREC LONG 0 ;rec number from start of file (cnt from 0)
TABLESIZE LONG 9130000 ;number of bytes of data in a record
DATAWIDTH INT 2 ;length of a single data sample in bytes
DATATYPE STRING '' ;type of data sample
NUMDIMS INT 2 ;number of dimensions in datatbl
DIM0 LONG 91300 ;number of data samples in 1st dimension
DIM1 LONG 100 ;number of ipps in 2nd dimension
NUMCHANNELS INT 1 ;number of chan. 1=single beam,2=dualbeam
TXSTART LONG 0 ;starting sample for tx in ipp
TXLEN LONG 11200 ;number of samples in tx window
DATASTART LONG 12000 ;starting sample in data window
DATALEN LONG 80000 ; number of samples in data window
NOISESTART LONG 92000 ;starting sample for noise (cal) window
NOISELEN LONG 100 ;number of samples in noise(cal) window
SYSTIME STRING '15:17:16' ;time for record (starttime ??)
AZ FLOAT 449.999;az position (deg) for rec (at start??)
ZAGREG FLOAT 8.46800;za (dome) position for rec
ZACH FLOAT 0.000200000;za (ch) position for rec.
;Notes:
; Above a Sample is an I or Q data sample. There are two data samples for
; each time sample.
; The startTime is the data sample if all of the samples were kept for
; output (they're not). There is usually a gap between the tx window
; and the data window and the datawindow and the cal window.
; To find the position in the data buffer for the data window you must
; skip over the txLen. To find the starting time for the data window you
; multiply the datastart *.5*sampletime (.5 since there are 2 data samples
; per time sample.
; shsget() seperates out the tx,data, and cal data for you.
;
; Input 10 records (of 100 ipps each)
;
--> istat=shsget(desc,d,nrec=10)
; in this case d is an array of 10 records:
; the structure contents are:
help,d,/
IDL> help,d,/st
DHDR STRUCT -> Array[1]; hods the data header
NIPPS LONG 100 ; number of ipps in record
D1 STRUCT -> Array[1]; the data
help,d.d1,/st
TX INT Array[2, 5600, 100] ;tx data:[i/q,5600txsmples ,100 ipps]
DAT INT Array[2, 40000, 100];hghtdat:[i/q,40000 smmples,100 ipps]
CAL INT Array[2, 50, 100] ;caldat :[i/q,50 calsmples ,100 ipps]
; --> position to a random record in the file
rec=35
istat=shspos(desc,rec)
NOTES:
- The shspos() routine assumes that all records are the same length
as the first record (this is supposed to be true).
(See /pkg/rsi/local/libao/phil/atm/0atmIntro.pro)
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)
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.
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(lun,spcBuf1,spcBuf2,nrec=nrec,baudLen=baudLen,
spclen=spclen,txSmpSkip=txSmpSkip,nheights=nheights,$
dinfo=dinfo,hdr=hdr,dotm=dotm,tmi=tmi,verbose=verbose
ARGS:
desc: struct desc that points at data file
KEYWORDS:
nrec: long number of records to process (default is 10).
Normally there are 100 ipps in a record so this would
be 10 seconds of data
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
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 nrec 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:
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).
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)
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 1).
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)
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)
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]
[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)