c RDR2TABLE is template for RDR formatted table output of chan. 1-5 c rdr2table reads the LOLA Reduced Data Record file and outputs c (0) spacecraft position, height, geoid, tx_pw, energy etc. c (1-5) longitude, latitude, height, pulsewidth, energy noise c **********************************************************************/ c revision history */ c v0.99 uses C routine unsigned.c: c V1.0 applies preliminary reflectance calibrations to the return energy/ c V1.1 dark energy subtracted / c V1.4 adds a boatload of options / c 2.0 version checks flags on bit 0 and bit 6 for consistency c 2.43 version implements W 5/05/2010 c 3.0 12/21/09 uses stream to read rdr-faster if whole file is output c 3.1 make oflag use consistent - starts time column at 0 if true c ************************************************************************ c double unsigned_ (unsigned int *s) return (double) *s; c gfortran -ffixed-line-length-none -O rdr2table.f unsigned.o iswap4.o kzext2.o spicelib.a -static-libgfortran -o rdr2table c ifort -132 -assume byterecl -O rdr2table.f unsigned.o iswap4.o kzext2.o -o rdr2table c for LITTLENDIAN version we swap words in lmet and tdt IMPLICIT NONE real*8 unsigned, d28, pi,r2d,d2r,r2de external unsigned parameter(d28=28.d0 + 1./8192.d0)! parameter(pi=3.141592653589794d0) parameter (r2d=180.d0/pi,d2r=pi/180.d0, r2de =9.d-03/pi) c ******* input record integer*4 rdr(64*7200*28), iwrd, mypos !stream integer*2 wdr(128*7200*28) equivalence(rdr,wdr) c ******* integer*8 LMET,TDT c these 8-byte integers equivalence(LMET,rdr(1))!LOLA MET Time as long long equivalence(TDT,rdr(3)) !Terrestrial dynamic time as long long integer*4 iswap4, itmp, numarg, ngrd c must be swapped integer*4 invalid /-2147483648/ integer*4 kzext2,zext external zext,kzext2 integer*4 MET,i,j,irec,orec,itdc,lseq,ifrm,iargc,subs,hz2fire real*8 offset,fract,time,geoid,xplse,xenrg,earle,earce,earpw real*8 solinc,solphs,offnadir,emission real*8 c32 /4294967296.d0/ real*8 et,sclkdp, orig /0./ real*8 sclon,sclat,scrad,ref real*8 radius/1737.4d0/ real*8 gain(5),thrs(5),bkgd(5),plse(5),rnge(5),enrg(5) real*8 dlon(5),dlat(5),drad(5) integer nois(5),gflg(5) real*8 rflct(5), trans(5),dark(5),fiber(5)!area-transm. product data dark /0.135d0,0.1159d0,0.1303d0,0.1775d0,0.2164d0/ integer*4 buf,rgstart,rgstop c minor frame counts as fractions of 5000000.d0 integer*4 minor(0:27) real*8 dscale/1.d-7/,gscale/1.d-6/, escale/1.d-6/ real*8 pscale/1.d-3/,rscale/1.d-6/, tscale/1.d-6/ character*16 arg character*28 utc character*64 filename! /'LOLARDR_00111N.DAT'/ logical swap /.false./,head/.false./,all/.false./,geop/.false./ logical flag/.false./,spot(0:5)/6*.false./,oflag/.false./ logical flagq/.false./,geographic/.true./,flagu/.false./,flagw/.false./ integer flagd /0/ !everything goes c fiber, and aft optics efficiency in percent data fiber /84.045654,87.437592,85.325625,90.273744,80.809245/ c diffractive optics splits energy into 5 spots, more in the center c these should be symmetrical but owing to measurement imprecision, c separate values are given in the cal rpt. This really needs to be flatfielded! data trans /0.22d0,0.14d0,0.14d0,0.14d0,0.14d0/!nominal c data trans /0.218d0,0.141d0,0.129d0,0.15d0,0.143d0/! c initialize optical link, transmission times rx telescope radius 0.07 squared do i=1,5 trans(i)=trans(i)*0.01d0*fiber(i)*0.0049d0 enddo c SPICE needs these if you uncomment et2utc c call ldpool('/Users/neumann/naif0009.tls') c initialize minor frame offsets - use to check range window against time minor(0)=0 do i=1,16 minor(i)=minor(i-1)+178571 enddo do i=17,27 minor(i)=minor(i-1)+178572 enddo numarg=iargc() !flag for headers, swapping if(numarg.eq.0) then write(*,*)'rdr2table 3.1-LOLA reduced data record to table file' write(*,*)'usage: rdr2table [0][A|1 2 3 4 5][C|c][D|d][F|f][G|g][[H|h](eader)][W|w][O|origin][[S|s](wap_bytes)]' goto 99999 endif do i=2,numarg call getarg(i,arg) if(arg(1:1).eq.'s' .or. arg(1:1).eq.'S') swap=.true. if(arg(1:1).eq.'h' .or. arg(1:1).eq.'H') head=.true. if(arg(1:1).eq.'f' .or. arg(1:1).eq.'F') flagd=255 ! any flag if(arg(1:1).eq.'m' .or. arg(1:1).eq.'M') flagd=255-64 !matlab flag only if(arg(1:1).eq.'n' .or. arg(1:1).eq.'N') flagd=254 !neumann's auto flag only c output bit 6 disagrees with bit 0 -failure of neumann's own? if(arg(1:1).eq.'q' .or. arg(1:1).eq.'Q') flagq=.true.!quarrel-matlab flag, auto flag c if(arg(1:1).eq.'u' .or. arg(1:1).eq.'U') flagu=.true.!not sure what this will mean if(arg(1:1).eq.'w' .or. arg(1:1).eq.'W') flagw=.true.!weird if(arg(1:1).eq.'c' .or. arg(1:1).eq.'C') flagd=flagd+2**10 !check for C-kernels if(arg(1:1).eq.'g' .or. arg(1:1).eq.'G') geop=.true. if(arg(1:1).eq.'o' .or. arg(1:1).eq.'O') oflag=.true. !MET origin if(arg(1:1).eq.'a' .or. arg(1:1).eq.'A') all =.true. if(arg(1:1).eq.'0' ) spot(0) =.true. if(arg(1:1).eq.'1'.or.all) spot(1) =.true. if(arg(1:1).eq.'2'.or.all) spot(2) =.true. if(arg(1:1).eq.'3'.or.all) spot(3) =.true. if(arg(1:1).eq.'4'.or.all) spot(4) =.true. if(arg(1:1).eq.'5'.or.all) spot(5) =.true. enddo c if (flagq.and.flagu) then c write(*,*) 'inconsistent usage of q(uarrel) and u(nflagged)' c stop c endif if(numarg.gt.0) call getarg(1,filename) c rdr? c open(unit=1,file=filename,recl=256,access='direct' c & ,status='old',err=9999) open(unit=1,file=filename,access="stream" & ,status='old',err=9999) read(1,end=3) rdr !get up to 7200 seconds 3 inquire(unit=1, POS=mypos)!points to next byte. close(1) irec=0 if(spot(0).and.head) & write(*,'(a,a)') &' SCLK_LOLA sc_alt 0 sc_longE sclat_N offnadir ' & ,'emission ifrm solinc solphs ngrd xenrg xplse' if(head.and.(spot(1).or.spot(2).or.spot(3).or.spot(4).or.spot(5))) & write(*,'(a,a)') &' SCLK_LOLA alt_km id longitudeE latitudeN range_km ' & ,' energy noise thrs gain flg reflect pulsewd' c if(swap) then MET=iswap4(rdr(1)) fract=unsigned(iswap4(rdr(2)))/c32 else MET=rdr(1) fract=unsigned(rdr(2))/c32 !cast fraction of met seconds to real*8 endif sclkdp = MET+fract c set origin time to MET, otherwise to start of orbit if(oflag) orig=sclkdp 1 continue j=0 iwrd=64*irec !offset to current record irec=irec+1 ! record number if(irec.ge.mypos/256) goto 999 c swap all words if(swap) then do i=1,60 rdr(iwrd+i)=iswap4(rdr(iwrd+i)) enddo c we should also swap halfwords at 61, 62, 63, 64 using kzext2 endif MET=rdr(iwrd+1) fract=unsigned(rdr(iwrd+2))/c32 !cast fraction of met seconds to real*8 c itmp=rdr(iwrd+4) !swap required for TDT c rdr(iwrd+4)=rdr(iwrd+3) c rdr(iwrd+3)=itmp geoid = rscale*rdr(iwrd+10) ref=radius if(geop) ref=geoid sclkdp = MET+fract c offset from 2000T12:00:00 to 2001T00:00:00 adopted by LRO is added c ET = tdt/c32 + 64.184d0 !plus a small periodic term c call et2utc (et,'isod',5,utc) c write(0,'(a,2f18.6,1x,a28)')'LMET,TDT',lmet/c32,tdt/c32,utc xenrg = escale*rdr(iwrd+5) xplse = pscale*rdr(iwrd+6) sclon = dscale*rdr(iwrd+7) if(sclon.lt.0.) sclon=sclon+360.d0 sclat = dscale*rdr(iwrd+8) scrad = rscale*unsigned(rdr(iwrd+9)) offnadir=r2de*kzext2(wdr(2*iwrd+121))!always found if(wdr(2*iwrd+122).ne. -1) then emission=r2de*kzext2(wdr(2*iwrd+122)) else emission=-1.d0! or zero? endif if(wdr(2*iwrd+123).ne. -1) then solinc=r2de*kzext2(wdr(2*iwrd+123)) else solinc=-1.d0 endif if(wdr(2*iwrd+124).ne. -1) then solphs=r2de*kzext2(wdr(2*iwrd+124)) else solphs=-1.d0 endif c truncate 28*fract ifrm = int(d28*fract) if(wdr(2*iwrd+128).ne. -1) then earpw=pscale*kzext2(wdr(2*iwrd+128)) else earpw=99.d0 endif ngrd=0 !if flagd mask is nonzero, this counts ground returns, otherwise any triggers do i=1,5 if(iand(rdr(iwrd+i*10+10),flagd).eq.0.and.rdr(iwrd+i*10+4).ne.-1) ngrd=ngrd+1 enddo c **********************************************************************/ if(spot(0))write(*,1005) sclkdp-orig,scrad-radius,0,sclon,sclat & ,offnadir,emission,ifrm,solinc,solphs,ngrd,xenrg,xplse c do i=1,5 flag =(iand(rdr(iwrd+i*10+10),flagd).eq.0) if(flagq)flag=flag.and.(iand(rdr(iwrd+i*10+10),64).ne.iand(rdr(iwrd+i*10+10),1)*64)!bit 6 differs from bit 0 if(flagw)flag=flag.or.(iand(rdr(iwrd+i*10+10),32768).ne.0)!bit 15 check output desired switches, valid geolocation, valid range if(flag.and.spot(i) & .and.rdr(iwrd+i*10+1).ne.invalid & .and.rdr(iwrd+i*10+4).ne.-1 & ) then dlon(i) = dscale*rdr(iwrd+i*10+1) if(geographic .and.dlon(i).lt.0.) dlon(i)=dlon(i)+360.d0 dlat(i) = dscale*rdr(iwrd+i*10+2) drad(i) = rscale*unsigned(rdr(iwrd+i*10+3)) rnge(i) = rscale*unsigned(rdr(iwrd+i*10+4)) plse(i) = pscale*rdr(iwrd+i*10+5) enrg(i) = escale*rdr(iwrd+i*10+6) nois(i) = rdr(iwrd+i*10+7) thrs(i) = tscale*rdr(iwrd+i*10+8) gain(i) = gscale*rdr(iwrd+i*10+9) gflg(i) = iand(255+1024+32768,rdr(iwrd+i*10+10)) c gflg(i) = iand( 1,rdr(iwrd+i*10+10)) c reflectivity if(flag .and. xenrg.gt.0.) then c trans(i) contains aperture, transmission efficiency, etc c scale factor accounts for km^2, mJ to fJ rflct(i)=rnge(i)**2 *escale*(enrg(i)-dark(i))/(xenrg*trans(i)) else rflct(i)=0.d0 endif write(*,1005) & sclkdp-orig,drad(i)-ref,i,dlon(i),dlat(i),rnge(i) & ,enrg(i),nois(i),thrs(i),gain(i),gflg(i),rflct(i) & ,plse(i) 1005 format & (f14.3,f10.4,i2,2f11.6,f10.3,f9.4,i6,f8.3,f8.3,i6,f10.5,f9.2) c endif enddo goto 1 999 continue close(1) goto 99999 9999 write(0,*)'file not opened: ',filename 99999 end