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 3.14 version implements L, adds 5 minutes to orbit, applies geographic
c 3.14159 fixes time limit test
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,offlo,offhi,zmax /180.d0/
	data offlo/0./,offhi/25.d0/
	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*7500*28), iwrd, mypos !stream
	integer*2 wdr(128*7500*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,lo,hi,lim
	real*8 dlonmin /-180.0/,dlonmax/360.d0/,dlatmin/-90.d0/,dlatmax/90.d0/,zenith
	real*8 timelo/0./,timehi/7500./,alt
	character*16 name(9)/'long_min','long_max','lat_min','lat_max','time_lo','time_hi','offnadlo','offnadhi','zenithmx'/,chr
	real*8 fred(9)
	equivalence(dlonmin,fred(1)),(dlonmax,fred(2)),(dlatmin,fred(3)),(dlatmax,fred(4)),(timelo,fred(5)),(timehi,fred(6)),(offlo,fred(7)),(offhi,fred(8)),(zmax,fred(9))
	real*8 offset,fract,time,geoid,xplse,xenrg,earle,earce,earpw
	real*8 solinc,solphs,offndr,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)
	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*128 filename
	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 SPICE needs these if you uncomment et2utc
c	if (flagu) call ldpool('/Users/neumann/naif0009.tls')
c initialize minor frame offsets - use to check range window against time
c	minor(0)=0
c	do i=1,16
c	 minor(i)=minor(i-1)+178571
c	enddo
c	do i=17,27
c	 minor(i)=minor(i-1)+178572
c	enddo
	numarg=iargc() !flag for headers, swapping
	if(numarg.eq.0) then
	 write(*,*)'rdr2table 3.14159-LOLA reduced data record to table file'
	 write(*,*)'usage: rdr2table <rdrfile>[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 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.!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. !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 ======= 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
	lo=timelo
	hi=timehi
	endif
	if(flagv) then
	do i=1,9
 	 write(*,6) name(i),' = ', fred(i)
	 write(*,*)
	enddo
	endif
6	 format(a,a,f10.3," ", $)
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
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=30) rdr  !get up to 7500 seconds
30	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
	lim = MET
	sclkdp = MET+fract
c set origin time to MET, otherwise to start of orbit
	if(oflag) orig=sclkdp
	ref=radius
c ================================================================== big loop
1       continue
	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
	else
c@@@@ smallendian TDT
	 itmp=rdr(iwrd+4)
	 rdr(iwrd+4)=rdr(iwrd+3)
	 rdr(iwrd+3)=itmp
c@@@@
	endif
	MET=rdr(iwrd+1)
c check time from start of orbit
	if (.not. limit .or.
     &   (MET - lim .ge.lo .and. MET -lim .lt. hi)
     &  ) then

	fract=unsigned(rdr(iwrd+2))/c32 !cast fraction of met seconds to real*8
	geoid = rscale*rdr(iwrd+10)
	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(geographic .and. sclon.lt.0.) sclon=sclon+360.d0
	sclat = dscale*rdr(iwrd+8)
	scrad = rscale*unsigned(rdr(iwrd+9))
	offndr=r2de*kzext2(wdr(2*iwrd+121))!always found
	zenith = r2de*kzext2(wdr(2*iwrd+123))
	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)) then
	if(.not.limit.or.
     &      (     sclon.ge.dlonmin .and. sclon.le. dlonmax
     &      .and. sclat.ge.dlatmin .and. sclat.le. dlatmax)
     &      .and.  offndr.ge.offlo   .and. offndr.le. offhi .and. zenith.le.zmax
     &  ) then

	write(*,1005) sclkdp-orig,scrad-radius,0,sclon,sclat
     & ,offndr,emission,ifrm,solinc,solphs,ngrd,xenrg,xplse
c
	endif !limit
	endif !spot 0
	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)

	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 )
     &  ) then

	 if(rdr(iwrd+i*10+3).ne.-1) then
	 drad(i) = rscale*unsigned(rdr(iwrd+i*10+3))
	  alt = max(-999.999d0,drad(i)-ref)
	 else
	  alt = -999.d0
	 endif
	 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,alt,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 !lon/lat limit
	 endif !flag
	enddo
	endif !met limits
        goto 1
999     continue
        close(1)
	goto 99999
9999	write(0,*)'file not opened: ',filename
99999   end
