;$Id: eis_data__readccsds.pro 405 2010-07-04 11:21:38Z viggoh $ ;+ ; NAME: ; EIS_DATA__READCCSDS ; ; PURPOSE: ; EIS_DATA__READCCSDS reads EIS data files and stores data and other ; information in the object fields of an EIS_DATA class object. ; ; CATEGORY: ; Hansteen/Wikstol Data analysis SW ; ; CALLING SEQUENCE: ; An object of class EIS_DATA is declared with this statement: ; eis_data = obj_new('eis_data') ; ; EIS_DATA__READ reads EIS data files and stores information ; into the object, with this statement: ; eis_data-> read, filename [, eis_hdr = eis_hdr] ; ; If called from within an eis_data procedure, the ; EIS_DATA__READ statement looks like this: ; self-> read, file [, eis_hdr = eis_hdr] ; ; INPUTS: ; file: Name of file (string) containing EIS data (and headers). ; ; KEYWORD PARAMETERS: ; eis_hdr: Set this keyword to obtain also the header ; information in the form of an object of type EIS_HDR. ; ; OUTPUTS: ; Reads EIS data files and fills the parameters of the EIS_DATA ; object with information. ; ; CALLS: ; EIS_HDR__DEFINE ; EIS_HDR__READ ; CCSDS_HDR__DEFINE ; CCSDS_HDR__READ ; LITTLE_ENDIAN ; ; COMMON BLOCKS: ; ; ; PROCEDURE: ; Description on the endian-tests: EIS CCSDS can come in the form of both ; little endian and big endian files. Depending on what endianess the ; computer that reads the data has, one may have to swap either header or ; pixel data or both. The fact that the header and pixel data are written ; to a CCSDS file as integers, while the header in EIS_HDR__READ is read ; as a byte array while the pixels are read as integers complicate ; the situation even more. Here is the "algorithm": ; IF Computer is big_endian: Swap both header and pixels if file is ; little endian. Swap nothing if file is big endian. ; IF computer is little endian: Swap header but not pixels if file is ; little endian. Swap pixels but not header if file is big endian. ; ; ; RESTRICTIONS: ; ; ; MODIFICATION HISTORY: ; 28-Mar-2001: Oivind Wikstol. ; 21-Sep-2001: Oivind Wikstol - Added documentation. ; 06-May-2004: Oivind Wikstol - self.hdr added ; 17-Jun-2004: Oivind Wikstol - Added self.filename ; 15-Feb-2005: Oivind Wikstol - Changed nraster added nexp_prp ; 26-Jan-2006: Oivind Wikstol - Fixed endian problem ; 16-Dec-2006: Viggo Hansteen - Added time of first and last ccsds ; packet ; 14-Jan-2007: Viggo Hansteen - Fixed bug introduced when attempting to; account for bad data packets that ; caused data with several exposures per ; raster position to be incorrectly read. ; 29-Nov-2007: A Gardini - Cleaned dangling references ; 11-Dec-2007: A Gardini - Reintroduction of hdr objects. ; 3-Jul-2008: A Gardini - Added init keyword. ; 18-Aug-2008: V Hansteen - Fixes to deal with corrupt or missing data ; 23-Sep-2009: V Hansteen - Further fixes to deal with ; missing or corrupt data ; 19-Feb-2010: V Hansteen - Yet further fixes to deal with ; corrupt or missing data, ; artificial fmirr and nexpc data added ; ; $Id: eis_data__readccsds.pro 405 2010-07-04 11:21:38Z viggoh $ ;- pro eis_data::readccsds, file, hdr = hdr, error = error, init=init if n_elements(init) eq 1 then return if n_params() ne 1 then begin message, 'No datafile specified!', /informational message, 'usage: read_data, file, hdr = hdr', $ /informational return endif self.title='Intensity' self.filename = file maxraster = 65000 ccsds_hdr_size = 10 eis_hdr_size = 256 ccsds_hdr = obj_new('ccsds_hdr') ; 10 bytes eis_buff = objarr(maxraster) ; 256 bytes data_len = lonarr(maxraster) hdrswap = 0 ; swap header bytes or not rstart: ; read start ; open the file counter = 0UL openr,lu, file, /get_lun fileinfo = fstat(lu) filesize = fileinfo.size ; get size of file in bytes ;check if filesize is zero. If so, print message and return ; if filesize eq 0 then begin ; print, 'File is empty!' ; return ; endif ; data = uintarr(filesize/2) ; int is two bytes. No more than this data in file ; ;(since header also takes up some bytes) i=0 repeat begin ; read first packet types (ccsds_hdr + eis_hdr) ccsds_hdr-> read, lu, swap = hdrswap ccsds_flag = ccsds_hdr-> getccsds_seq_flag() self.ccsds_packet_time.last=ccsds_hdr->getccsds_time() ; 8 bytes added to data in Japan in october 2004. Set ver_no to ; determine how to read data pkt_len = ccsds_hdr-> getccsds_pkt_len() if (pkt_len lt 260) then ver_no = 0 else ver_no = 1 eis_buff[i] = obj_new('eis_hdr') eis_buff[i]-> setver_no, ver_no ; read md hdr eis_buff[i]-> read, lu, swap = hdrswap ; first time - check if file is big-endian or little endian: ; if little endian: set swap keyword to 1 and start all over if i eq 0 then begin self.ccsds_packet_time.first=ccsds_hdr->getccsds_time() datatype = (eis_buff[i]-> getdatatype()) if datatype ne 195 then begin if little_endian() then begin hdrswap = 1 close,lu & free_lun,lu obj_destroy,eis_buff[0] goto, rstart endif else begin hdrswap = 1 close,lu & free_lun,lu obj_destroy,eis_buff[0] goto, rstart endelse endif else begin if little_endian() then begin if (hdrswap) then pixswap = 0 else pixswap = 1 endif else begin pixswap = hdrswap endelse endelse endif ; compute data size if i eq 0 then begin exp_info=eis_buff[i]->getexp_info() nwin=exp_info.nwin ; number of line windows nexp=exp_info.nexp ; number of exposures nexpc=exp_info.nexpc nexp_prp = exp_info.exp_prp fmirr_step=exp_info.fmir_step fmirr0=exp_info.fmirr_pos-((nexpc-1)/nexp_prp)*fmirr_step fmirr=(findgen(nexp/nexp_prp)*fmirr_step+fmirr0)#(intarr(nexp_prp)+1) fmirr=fmirr[sort(fmirr)] ; print,fmirr yw=exp_info.yw ; window height npix=0L ; size of exposure for iw=0,nwin-1 do npix=long(yw)*long(exp_info.x_l[iw])+npix datasize=npix*nexp data=uintarr(datasize) endif ; read data type packets (ccsds_hdr + pixel data) if not eof(lu) then begin ; read all data packets for 1 raster position repeat begin ccsds_hdr-> read, lu, swap = hdrswap self.ccsds_packet_time.last=ccsds_hdr->getccsds_time() ; subtract sizeof(sec. hdr) (4 bytes) packet_len = (ccsds_hdr-> getccsds_pkt_len()-4+1)/2 pxd=uintarr(packet_len) data_len[i] = long(packet_len) + data_len[i] ; read pixel data into the 'buf' variable readu, lu, pxd if pixswap ne 0 then byteorder, pxd ; pixel data is 14 bits: pxd=ishft(ishft(pxd,2),-2) ; store buf into the data array: data(counter:counter+packet_len-1) = pxd counter = counter+packet_len ccsds_flag = ccsds_hdr-> getccsds_seq_flag() ; kluge to account for "bad" data set (was "ne 0") ; normally flag should be 1 for first packets ; 0 for continuation packets and ; 2 for last packets ; ; in md_data.3 and md_data.4 flags seem to be incremented by 1. endrep until ccsds_flag ge 2 endif i=i+1 endrep until eof(lu) ; store all i eis_hdrs ; ; check how many packets are missing and which exposures are complete ; nexp=(eis_buff[0]->getexp_info()).nexp message,'assuming nexp = '+string(nexp,format='(i4)')+' exposures in file',/info exp_ok=intarr(nexp) isub=0 size_packets=lonarr(eis_buff[0]->getnr_of_packets()) nexpc=(eis_buff[0]->getexp_info()).nexpc for ip=0,i-1 do begin jsub=eis_buff[ip]->getsub_seq_count() nr_pack=eis_buff[ip]->getnr_of_packets() if nexpc ne (eis_buff[ip]->getexp_info()).nexpc then begin nexpc=(eis_buff[ip]->getexp_info()).nexpc isub=0 endif ; print,ip,isub,jsub,nr_pack,(jsub+1) mod nr_pack,nexpc,eis_buff[ip]->getserial_packet_no() if ( (jsub+1) mod nr_pack eq 0 ) and ( isub eq jsub) then begin exp_ok[nexpc-1]=1 isub=-1 endif isub=isub+1 endfor ip=0 ok=0 offset=long((eis_buff[ip]->getexp_info()).nexpc-1)*(eis_buff[ip]->getnr_of_packets()) repeat begin if exp_ok[(eis_buff[ip]->getexp_info()).nexpc-1] then begin first_good=eis_buff[ip]->getserial_packet_no() for isub=0,eis_buff[ip]->getnr_of_packets()-1 do size_packets[isub]=data_len[ip+isub] ok=1 offset=((eis_buff[ip]->getexp_info()).nexpc-1)*(eis_buff[ip]->getnr_of_packets()) ; endif else offset=offset+1L endif ip=ip+1 endrep until ok or ip gt i-1 if (not ok) and (ip gt i-1) then begin message,'no complete exposure found in data!!!',/info message,'computing offset to missing packet based on sub sequence count',/info message,'data may be offset!!!',/info first_good=eis_buff[0]->getserial_packet_no() offset=long((eis_buff[0]->getexp_info()).nexpc-1)*(eis_buff[0]->getnr_of_packets()) offset=eis_buff[0]->getsub_seq_count()+offset endif first=first_good-offset last=first+eis_buff[0]->getnr_of_packets()*nexp-1 if (where(exp_ok eq 0))[0] ne -1 then $ message,'exposure(s) '+strjoin(string(where(exp_ok eq 0)+1,format='(i4)'))+' are missing or damaged',/info nhdr=last-first+1 hdr=objarr(nhdr) for j=0,i-1 do begin hdr[eis_buff[j]->getserial_packet_no()-first]=eis_buff[j] endfor pointer=0L for j=0,nhdr-1 do begin size_packet=size_packets[j mod eis_buff[0]->getnr_of_packets()] if not obj_valid(hdr[j]) then begin hdr[j]=obj_new('eis_hdr') serial_packet_no=first+j exp_info.nexpc=j/(eis_buff[0]->getnr_of_packets())+1 exp_info.ti_1=0L exp_info.ti_2=0L exp_info.exp_mhc=0UL exp_info.exp_dur=0U exp_info.fmirr_pos=fmirr[exp_info.nexpc-1] hdr[j]->setexp_info,exp_info hdr[j]->setcorrupt,1 hdr[j]->setmain_id,eis_buff[0]->getmain_id() hdr[j]->setserial_packet_no,serial_packet_no hdr[j]->setnr_of_packets,eis_buff[0]->getnr_of_packets() hdr[j]->setbit_compr_mode,eis_buff[0]->getbit_compr_mode() hdr[j]->setimg_compr_mode,eis_buff[0]->getimg_compr_mode() hdr[j]->sethuffman_ac,eis_buff[0]->gethuffman_ac() hdr[j]->sethuffman_dc,eis_buff[0]->gethuffman_dc() hdr[j]->setquant_tab,eis_buff[0]->getquant_tab() hdr[j]->setline_id,eis_buff[0]->getline_id() data[pointer:*]=shift(data[pointer:*],size_packet) endif ; print,(hdr[j]->getexp_info()).fmirr_pos,(hdr[j]->getexp_info()).nexpc pointer=pointer+size_packet endfor ;print,data ; hdr is also part of data object self.hdr=ptr_new(obj_new()) *self.hdr = hdr[offset] ; obj_destroy,eis_buff obj_destroy,ccsds_hdr ; fill each line window with data, and set xstart, xlength and ystart ; for each line window exp_info=hdr[offset]->getexp_info() nwin=exp_info.nwin ; number of line windows self.nexp = (hdr[offset]-> getexp_info()).nexp ; number of exposures ; number of exposures pr raster position: self.nexp_prp = (hdr[offset]-> getexp_info()).exp_prp self.nraster = self.nexp/self.nexp_prp self.nwin=nwin self.ys=exp_info.yws self.yw=exp_info.yw ; distinguish between sit-and-stare and scanning: fmirr_ss = (hdr[offset]-> getexp_info()).fmir_step if fmirr_ss eq 0 then self.obsmode = 'sit-and-stare' else $ self.obsmode = 'scanning' ; determine if sit-and-stare or rastering: if exp_info.fmir_step eq 0 then begin self.sit_and_stare = 1 self.nraster = 1 endif for i=0,nwin-1 do begin addpix=0 ; according to the new definitions of the CCD's, the long wavelength CCD ; is now called A, and the short wavelength CCD is called B. Therefore, the ; CCD B (short) is associated with CCD flags 2 and 3 (10 and 11) and CCD A (long) ; has CCD flags 0 and 1 (00 01). This means we add 2148 pixels if ccd_nr is lt 2. ; OW, 5/12/01, after talking to Khalid and Ady. ; ;;; if exp_info.ccd_nr[i] gt 1 then addpix=2148 ; earlier statement (OW, 5/12/01) ; if exp_info.ccd_nr[i] lt 2 then addpix=2148 self.xs[i]=exp_info.x_s[i]+addpix self.xw[i]=exp_info.x_l[i] self.w[i]=ptr_new(uintarr(self.xw[i],self.yw[i],self.nexp)) endfor self.exp = ptr_new(fltarr(self.nexp)) self.aux_data.fmirr = ptr_new(uintarr(self.nexp)) self.aux_data.mhc_dur = ptr_new(ulonarr(self.nexp)) nrppexp = hdr[offset]->getnr_of_packets() ; number of packets pr. exposure for j=0,self.nexp-1 do begin (*self.exp)[j] = (hdr[j*nrppexp]->getexp_info()).exp_dur/1.e2 ; exposure time in seconds (*self.aux_data.fmirr)[j] = (hdr[j*nrppexp]->getexp_info()).fmirr_pos (*self.aux_data.mhc_dur)[j] = (hdr[j*nrppexp]->getexp_info()).exp_mhc endfor jumps = lonarr(nwin, self.nexp) rastersz = 0l for i = 0, nwin-1 do begin rastersz = long(self.xw[i])*long(self.yw[i])+rastersz endfor for j = 0, self.nexp-1 do begin jumps[0, j] = rastersz*j for i = 1, nwin-1 do begin jumps[i, j] = jumps[i-1, j]+long(self.xw[i-1])*long(self.yw[i-1]) endfor endfor last_pix = n_elements(data) for i=0,nwin-1 do begin xw = long(self.xw[i]) yw = long(self.yw[i]) rd=uintarr(xw,yw,self.nexp) for j=0,self.nexp-1 do begin ; if (hdr[j*nrppexp]->getexp_info()).ti_1 ne 0 then begin ; pknr=hdr[j*nrppexp]->getserial_packet_no() ; jexp=(pknr-first_good)/nrppexp ; print,jexp ; dl = n_elements(data) ; if jumps[i,j]+xw*yw gt dl then begin ; rd[*,*,j] = 0U ; endif else begin wd=reform(data[jumps[i, j]:jumps[i, j]+xw*yw-1],yw, xw) rd[*,*,j]=transpose(wd) ; endelse ; endif endfor *self.w[i]=rd endfor close,lu & free_lun,lu end