c RDR2TAB 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 Q checks flags on bit 0 and bit 6 for consistency c 2.42 version implements Q 4/15/2010 c 2.43 version implements W 5/05/2010 c 2.44 version implements L 5/30/2010 c 2.45 version implements D 7/30/2010 c 2.46 version implements U 9/30/2010 c renamed rdr2tab 8/30/2010 c ************************************************************************ c double unsigned_ (unsigned int *s) return (double) *s; c gfortran -ffixed-line-length-none -O rdr2tab.f unsigned.o iswap4.o kzext2.o -static-libgfortran spicelib.a -o rdr2tab c for LITTLENDIAN version we swap words in lmet and tdt IMPLICIT NONE real*8 unsigned, d28, pi,r2d,d2r,r2de,utcsec,del2 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) integer*2 wdr(128) equivalence(rdr,wdr) c ******* integer*8 LMET,TDT c these 8-byte integers equivalence(LMET,rdr(1))!not currently used equivalence(TDT,rdr(3)) 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,lo,hi,lim real*8 dlonmin /-180.0/,dlonmax/360.d0/,dlatmin/-90.d0/,dlatmax/90.d0/ real*8 zenith, zmax/180.d0/,zmin/0./ real*8 timelo/0./,timehi/7200./ real*8 offndr, offlo/0./,offhi/21.d0/ real*8 offset,fract,time,geoid,xplse,xenrg,earle,earce,earpw real*8 solinc,solphs,emission real*8 c32 /4294967296.d0/ real*8 tmp,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),longflag 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*22 utc character*128 filename! /'LOLARDR_093330622.DAT'/ character*128 limitfile /'r2x.prm'/ logical swap /.false./,head/.false./,all/.false./,geop/.false./ logical flag/.false./,spot(0:5)/6*.false./,oflag/.false./,limit/.false./ logical flagq/.false./,geographic/.true./,flagu/.false./,flagv/.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 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(*,*)'rdr2tab 2.46 - LOLA reduced data record to table file NOT using stream i/o' write(*,*)'usage: rdr2tab [A 0 1 2 3 4 5][F|flag][G|geoid][H|header][L|limitfile][O|origin]' write(*,*)'[C|ckernel][D|dateline][M|manual][N|n][O|origin][Q|quarrel][U|utc][[S|swap_bytes][V|verbose][W|weird]' 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.'d' .or. arg(1:1).eq.'D') geographic=.false. !dateline 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.'l' .or. arg(1:1).eq.'L') then limit =.true. if(len_trim(arg).gt.1) limitfile = arg(2:) endif if(arg(1:1).eq.'m' .or. arg(1:1).eq.'M') flagd=255-64 !matlab flag overrides auto 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 if(arg(1:1).eq.'u' .or. arg(1:1).eq.'U') flagu=.true.!utc output if(arg(1:1).eq.'v' .or. arg(1:1).eq.'V') flagv=.true. !verbose 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. 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 SPICE needs these if you uncomment et2utc if (flagu) call ldpool('naif0009.tls') c ======= limits if(limit) then open(unit=1,file=limitfile,status='old',err=3) read(1,*,end=2) dlonmin,dlonmax,dlatmin,dlatmax,timelo,timehi,offlo,offhi,zmax!,zmin 2 close(1) goto 4 3 continue write (0,*)'limit file [e.g., r2x.prm] not found' stop 4 continue if(dlonmin.ge.dlonmax) then tmp=dlonmax dlonmax=dlonmin dlonmin=tmp endif if(dlatmin.ge.dlatmax) then tmp=dlatmax dlatmax=dlatmin dlatmin=tmp endif if(zmin.ge.zmax) then tmp=zmax zmax=zmin zmin=tmp endif lo=timelo hi=timehi 6 format(a,a,f10.3," ", $) endif ref=radius check if output specified if(.not.(spot(0).or.spot(1).or.spot(2).or.spot(3).or.spot(4).or.spot(5)))then write(0,*)'Must use at least one option (0, 1, 2, 3, 4, 5, or a)' goto 99999 endif if(numarg.gt.0) call getarg(1,filename) c rdr? open(unit=1,file=filename,recl=256,access='direct' & ,status='old',err=9999) irec=0 if(flagv) write(*,*)'flagd',flagd 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 ================================================================== big loop 1 continue j=0 irec=irec+1 read(1,rec=irec,err=999) rdr c swap all words if(swap) then do i=1,60 rdr(i)=iswap4(rdr(i)) enddo else c@@@@ smallendian itmp=rdr(4) rdr(4)=rdr(3) rdr(3)=itmp c@@@@ endif c ===== time MET=rdr(1) if (irec.eq.1) lim = MET fract=unsigned(rdr(2))/c32 !cast fraction of met seconds to real*8 if(flagu) then c offset from 2000T12:00:00 to 2001T00:00:00 adopted by LRO is already added to MET c utcsec =(tdt-tai2tdt-leap)/c32 in 8-byte arithmetic utcsec = tdt/c32-66.184d0 ! difference between TDT and UTC call DELTET (utcsec, 'UTC', DEL2)! Return the value of Delta ET (ET-UTC) c write(0,*)'del2',del2 et = utcsec + DEL2 call et2utc (et,'isod',3,utc) c call et2utc (utcsec + DEL2,'isod',3,utc) endif ! utc time geoid = rscale*rdr(10) if(geop) ref=geoid c itmp=rdr(4) !tdt long integer is swapped c rdr(4)=rdr(3) c rdr(3)=itmp 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 sclkdp = MET+fract !seconds, not ticks of course if(oflag.and.irec.eq.1) orig=sclkdp !start of pass c check time from start of orbit if (.not. limit .or. & (MET - lim .gt.lo .and. MET -lim .lt. hi) & ) then c xenrg = escale*rdr(5) xplse = pscale*rdr(6) sclon = dscale*rdr(7) if(sclon.lt.0.) sclon=sclon+360.d0 sclat = dscale*rdr(8) scrad = rscale*unsigned(rdr(9)) offndr=r2de*kzext2(wdr(121))!always found zenith = r2de*kzext2(wdr(123)) if(wdr(122).ne. -1) then emission=r2de*kzext2(wdr(122)) else emission=-1.d0! or zero? endif if(wdr(123).ne. -1) then solinc=r2de*kzext2(wdr(123)) else solinc=-1.d0 endif if(wdr(124).ne. -1) then solphs=r2de*kzext2(wdr(124)) else solphs=-1.d0 endif c truncate 28*fract ifrm = int(d28*fract) c earle=dscale*unsigned(rdr(63)) !V.2.2 c earce=unsigned(rdr(63))/c32!fraction of second if(wdr(128).ne. -1) then earpw=pscale*kzext2(wdr(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(i*10+10),flagd).eq.0.and.rdr(i*10+4).ne.-1) ngrd=ngrd+1 enddo c **********************************************************************/ if (.not.flagu) write(unit=utc,fmt='(f22.6)') sclkdp-orig if(spot(0))write(*,1005) utc(6:), scrad-radius,0,sclon,sclat & ,offndr,emission,ifrm,solinc,solphs,ngrd,xenrg,xplse c do i=1,5 longflag=rdr(i*10+10) !long flag word c check output for ground, geolocation, valid range flag =(iand(longflag,flagd).eq.0).and.(rdr(i*10+1).ne.invalid).and.(rdr(i*10+4).ne.-1) if(flagq)flag=flag.and.(iand(longflag,64).ne.iand(longflag,1)*64)!bit 6 differs from bit 0 if(flagw) flag=iand(longflag,flagd).eq.0.or.(iand(longflag,32768).ne.0)!bit 15 check output desired switches, valid geolocation, valid range if(flag.and.spot(i) & ) then dlon(i) = dscale*rdr(i*10+1) if(geographic .and.dlon(i).lt.0.) dlon(i)=dlon(i)+360.d0 dlat(i) = dscale*rdr(i*10+2) c if(flagv.and.limit)write(*,*)i,dlonmin,dlonmax,dlatmin,dlatmax,offlo,offhi,zmax,zmin if(.not.limit.or. & ( dlon(i).ge.dlonmin .and. dlon(i).le. dlonmax & .and. dlat(i).ge.dlatmin .and. dlat(i).le. dlatmax & .and. offndr.ge.offlo .and. offndr.le. offhi & .and. zenith.le.zmax .and. zenith .ge.zmin ) & ) then drad(i) = rscale*unsigned(rdr(i*10+3)) rnge(i) = rscale*unsigned(rdr(i*10+4)) plse(i) = pscale*rdr(i*10+5) enrg(i) = escale*rdr(i*10+6) nois(i) = rdr(i*10+7) thrs(i) = tscale*rdr(i*10+8) gain(i) = gscale*rdr(i*10+9) c gflg(i) = rdr(i*10+10) c bits 0-7,10,15 gflg(i) = iand(255+1024+32768,rdr(i*10+10)) c gflg(i) = iand( 1,rdr(i*10+10)) c reflectivity if(flag .and. xenrg.ge. 1.d0 ) then ! should be 2.5 to 3.2 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) & utc(6:),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 & (a,f10.4,i2,2f11.6,f10.3,f9.4,i6,f8.3,f8.3,i6,f10.5,f9.2) c endif endif !dlon/lat limits enddo !1-5 endif !met limits goto 1 999 continue close(1) goto 99999 9999 write(0,*)'file not opened: ',filename 99999 end