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 2.5 fixed time limit test
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/7500./,alt
	real*8 offndr, offlo/0./,offhi/25.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.5 - LOLA reduced data record to table file NOT using stream i/o'
	 write(*,*)'usage: rdr2tab <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 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
	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
        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 .ge.lo .and. MET -lim .lt. hi)
     &  ) then
c
	xenrg = escale*rdr(5) 
	xplse = pscale*rdr(6) 
	sclon = dscale*rdr(7)
	if(geographic .and. 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))then
	if(.not.limit.or.
     &      (     sclon.ge.dlonmin .and. sclon.le. dlonmax
     &      .and. sclat.ge.dlatmin .and. sclat.le. dlatmax )
     &  ) then
	write(*,1005) utc(6:), scrad-radius,0,sclon,sclat
     & ,offndr,emission,ifrm,solinc,solphs,ngrd,xenrg,xplse
	endif !limit
	endif !spot 0
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)
	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

	 if(rdr(i*10+3).ne.-1) then
	 drad(i) = rscale*unsigned(rdr(i*10+3))
	  alt = max(-999.999d0,drad(i)-ref)
	 else
	  alt = 0.
	 endif
	 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:),alt,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
