c RDR2XYZ, Greg Neumann, GSFC
c          edited by Kopal Jha, Sigma/GSFC
c
c Read in an RDR, and write out a binary file (or set of files)
c corresponding to the channels (0-5) requested.  The data can be limited by 
c parameters such as lat and long, so that data is written out only when 
c appropriate.
c
c rdr2xyz 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 ************************************************************************
c
c revision history starting with rdr2table */
c
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 for albedo /
c V2.0 file output argument required terminating at "." /
c V2.4 option "l" checks for offnadir angle limits as well as lat/lon
c 2.42 option "f" doesn't check ground, used for spot 1 to get ground track
c 2.43 version implements W 5/05/2010 to allow 'weird' points
c 3.0 uses stream to read rdr-faster
c 3.1 adds solar zenith
c 3.14 uses open status = 'replace' for output, close status = 'delete' if empty
c 3.14159 parses option "l" for limitfile
c 3.1416 increase array size, fix time limit test
c 3.what? xmit enrgy, not pulsewidth
c 4.6692 s/c height
c*******************************************************************************
c
c double unsigned_ (unsigned int *s) return (double) *s;
c
c To complile:
c gfortran -ffixed-line-length-none -static-libgfortran -O -o rdr2xyz rdr2xyz.f 
c          unsigned.o iswap4.o kzext2.o
c
c          .o files here compiled using gcc -cO
c
c optionally  spicelib.a
c for LITTLENDIAN version we swap words in lmet and tdt
c
c*******************************************************************************

	IMPLICIT NONE
	real*8 unsigned, d28, pi,r2d,r2de,offndr/0./,offlo,offhi,zenith, data(5),height
	external unsigned
	parameter(d28=28.d0 + 1./8192.d0)!
	parameter(pi=3.141592653589794d0)! cancels out for albedo
	parameter (r2d=180.d0/pi,r2de=9.d-03/pi) !radians to degrees, short int to degrees
	integer*4 rdr(64*7500*28), iwrd, mypos !stream
	integer*2 wdr(128*7500*28)
	equivalence(rdr,wdr)
	character*6 disp /'keep  '/, version /'4.6692'/
	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))
	integer*4  lim, lo, hi
	data lim /0/,lo/0/,hi/7500/, offlo/0./,offhi/25.d0/
	integer*8 LMET,TDT
c these 8-byte integers must be swapped
c	equivalence(LMET,rdr(iwrd+1))!not currently used
c	equivalence(TDT,rdr(iwrd+3))
c if they are used
	integer*4 iswap4, itmp, numarg, kzext2
	integer*4 invalid /-2147483648/
	integer*4 jzext,zext
	external zext,jzext
	integer*4 irec1(5),idot
	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 c32 /4294967296.d0/
	real*8 et,sclkdp,orig /0./
	real*8 sclon,sclat,scrad,ref,tmp
	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 albed(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/
	real*8 dlonmin /-180.0/,dlonmax/360.d0/,dlatmin/-90.d0/,dlatmax/90.d0/
	real*8 zmax/180.d0/,zmin/0./,timelo/0./,timehi/7500./
	character*16 arg
	character*28 utc
        character*128 filename! /'LOLARDR_09180000.DAT'/
        character*128 fileout ! /'LOLARDR_09180000.xyz'/
	character*128 limitfile /'r2x.prm'/
	logical swap/.false./,bidf/.false./,flagp/.false./,geop/.false./,limit/.false./
	logical flag/.false./,hflag/.false./,rflag/.false./,xflag/.false./,yflag/.false./,all/.false./,tflag/.false./
	integer flagd /255/
	logical spot(0:5)/6*.false./,extend/.false./,flagq/.false./,geographic/.true./,flagv/.false./,flagw/.false./
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
	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/!14.1,12.9,15.0 and 14.3
c
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
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
c                  
c   get command line arguments
c
	numarg=iargc() !flag for headers, swapping
        call getarg(1,filename)
        call getarg(2,fileout)
	if(numarg.lt.2) then
	 write(*,'(a,a,a)') 'rdr2xyz ',version,'-LOLA reduced data record to binary files'
	 write(*,*)'usage: rdr2xyz <rdrfile outfileroot [options]>'
	 write(*,*)'  [1-5] selected spots'
	 write(*,*)'  Invalid ranges are excluded by default (unless F)'
	 write(*,*)'  [A|a]ll 5 spots'
	 write(*,*)'  [B|b]i-directional reflectance'
	 write(*,*)'  [D|d]ateline (plus/minus 180 instead of 0-360)'
	 write(*,*)'  [E|e]xtend beyond 0, 360 (if lon+360 < dlonmax)'
	 write(*,*)'  [F|f]lag for ground is NOT checked'
	 write(*,*)'  [G|g]eoid reference instead of IAU radius (1737.4 km)'
	 write(*,*)'  [H|h]eight (above or below reference datum)-Default'
	 write(*,*)'  [L|l]imit_file <r2x.prm> (lonmin, lonmax, latmin, latmax, tmin, tmax, offndr, zmax)'
	 write(*,*)'  [M|m]atlab manual flag bit only test. Default is M and N'
	 write(*,*)'  [N|n]eumanns own flag bit - automatic density, slope test'
	 write(*,*)'  [P|p]ulsewidth in nanoseconds (full width)'
	 write(*,*)'  [Q|q]uarrel: mat flag, auto flag disagree'
	 write(*,*)'  [V|v]erbose: debug output'
	 write(*,*)'  [R|r]ange to surface(slant range)'
	 write(*,*)'  [S|s]wap bytes if bigendian data on PC or vice versa'
	 write(*,*)'  [T|t]riggers, output if any flag bits are SET'
	 write(*,*)'  [W|w]eird flags are considered ground'
	 write(*,*)'  [X|x]transmit energy (mJ) on 0'
	 write(*,*)'  [Y|y]transmit pulsewidth (ns) on 0'
	 goto 99999
	endif
c                  
c                  read in command line arguments
c
	do i=3,numarg
	 call getarg(i,arg)
	 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.'s' .or. arg(1:1).eq.'S') swap=.true.
	 if(arg(1:1).eq.'g' .or. arg(1:1).eq.'G') geop=.true.
	 if(arg(1:1).eq.'a' .or. arg(1:1).eq.'A') all =.true.
	 if(arg(1:1).eq.'b' .or. arg(1:1).eq.'B') bidf=.true. !BidirReflect
	 if(arg(1:1).eq.'d' .or. arg(1:1).eq.'D') geographic=.false. !non-dateline
	 if(arg(1:1).eq.'e' .or. arg(1:1).eq.'E') extend=.true.!
	 if(arg(1:1).eq.'f' .or. arg(1:1).eq.'F') flagd=255 !default
	 if(arg(1:1).eq.'h' .or. arg(1:1).eq.'H') hflag=.true. !height-default
	 if(arg(1:1).eq.'m' .or. arg(1:1).eq.'M') flagd=255-64 !manual flag only
	 if(arg(1:1).eq.'n' .or. arg(1:1).eq.'N') flagd=254 !neumann's flag only
	 if(arg(1:1).eq.'p' .or. arg(1:1).eq.'P') flagp=.true. !pulsewidth
	 if(arg(1:1).eq.'q' .or. arg(1:1).eq.'Q') flagq=.true. !quarrel-manual flag disagrees with auto flag
	 if(arg(1:1).eq.'v' .or. arg(1:1).eq.'V') flagv=.true. !verbose
	 if(arg(1:1).eq.'r' .or. arg(1:1).eq.'R') rflag=.true. !range
	 if(arg(1:1).eq.'t' .or. arg(1:1).eq.'T') tflag=.true. !test for bad ground flags
	 if(arg(1:1).eq.'x' .or. arg(1:1).eq.'X') xflag=.true. !transmit energy on 0
	 if(arg(1:1).eq.'y' .or. arg(1:1).eq.'Y') yflag=.true. !transmit pulsewd on 0
	 if(arg(1:1).eq.'w' .or. arg(1:1).eq.'W') flagw=.true.!weird
	 if(arg(1:1).eq.'0'.or.xflag.or.yflag) 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

	idot=index(fileout,".")
	if (idot.le.1) then
	  write(*,*) 'outfileroot must contain "."'
	  stop 42
	endif

c                       open named RDR file, skip out if does not exist
c open rdr?

	open(unit=1,file=filename,access="stream"
     & ,status='old',err=9999)
	read(1,end=1) rdr  !get up to 7200 seconds
1	inquire(unit=1, POS=mypos)!points to next byte.
	close(1)

        irec=0
	orec=0

	if(spot(0))open(unit=10,file=fileout(1:idot)//'0.xyz',access='direct',recl=24,status='replace')
	if(spot(1))open(unit=11,file=fileout(1:idot)//'1.xyz',access='direct',recl=24,status='replace')
	if(spot(2))open(unit=12,file=fileout(1:idot)//'2.xyz',access='direct',recl=24,status='replace')
	if(spot(3))open(unit=13,file=fileout(1:idot)//'3.xyz',access='direct',recl=24,status='replace')
	if(spot(4))open(unit=14,file=fileout(1:idot)//'4.xyz',access='direct',recl=24,status='replace')
	if(spot(5))open(unit=15,file=fileout(1:idot)//'5.xyz',access='direct',recl=24,status='replace')
c
c    records for 5 channels
	do i=1,5
	 irec1(i)=0
	enddo
	if(flagv) write(*,*)'fileout = ',fileout(1:idot)//'?.xyz'
c default is to read limits from r2x.prm, if not, stdinput.
c for anything beyond a few rdrs, this parameter file becomes essential

	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)
	if(flagv) then
	do i=1,9
 	 write(*,6) name(i),' = ', fred(i)
	 write(*,*)
	enddo
	endif
	goto 4
3	continue
        write (*,*)'limit file [e.g., r2x.prm] not found, waiting for keyboard input'
	do i=1,9
 	 write(*,6) name(i),' = ', fred(i)
	 read(*,'(a)',end=301) chr
	 read(unit=chr,fmt=*,end=301) fred(i)
301	 continue
	enddo
	open(unit=1,file='r2x.prm',status='new',err=4)
	do i=1,9
	 write(1,*) fred(i)
	enddo
	close(1)
4	continue
	lo = timelo
	hi = timehi
	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
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
c **********************BIGLOOP*******************************************/
100     continue
	iwrd=64*irec !offset to current record
        irec=irec+1 ! record number
	if(irec.ge.mypos/256) goto 999
c swap all words - machine specific, and only if flagged
	if(swap) then
	 do i=1,60
	  rdr(iwrd+i)=iswap4(rdr(iwrd+i))
	 enddo
c halfwords at 61 - 64 are accessed with function kzext2, which should perform the swapping
	endif
c
	MET=rdr(iwrd+1)
	if (irec.eq.1) lim = MET
	if (.not. limit .or.
     &   (MET - lim .ge.lo .and. MET -lim .lt. hi)
     &  ) then
	offndr = r2de*kzext2(wdr(2*iwrd+121))
	zenith = r2de*kzext2(wdr(2*iwrd+123))
	if(geop) ref=rscale*rdr(iwrd+10) !geoid
	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))
c truncate 28*fract
c 	ifrm = int(d28*fract)
c	earle=dscale*unsigned(rdr(iwrd+63))
c	earce=unsigned(rdr(iwrd+63))/c32!fraction of second
c	earpw=pscale*kzext2(wdr(2*iwrd+128))
c **********************************************************************/
	if(xflag) then
	 orec=orec+1
	 write(10,rec=orec)sclon,sclat,xenrg
	else if (yflag) then
	 orec=orec+1
	 write(10,rec=orec)sclon,sclat,xplse
	else if (spot(0)) then !default
	 orec=orec+1
	 height = rscale*unsigned(rdr(iwrd+9)) - ref
	 write(10,rec=orec)sclon,sclat,height
	endif
c output spot loop -test for ground shots or shots matching flag bits
	do i=1,5
c eliminated
c     & .and.rdr(iwrd+i*10+3).ne.-1 !radius may be invalid, but spot 1 is still geolocated 
c	 if(flagq) flag = bit 6 clear, bit 0 set
c	 if(flagu) flag = bit 6 set, bit 0 clear?
cc	 if(flagq.and.(iand(rdr(iwrd+i*10+10),64).ne.iand(rdr(iwrd+i*10+10),1)*64)) then !bit 6 differs from bit 0
cq	   flagd=0 !all good below
cq	 else
cq	   rdr(iwrd+i*10+10) =-1!reject
cq	 endif
cq	 if(flagv) write(*,*)'flagd',flagd
	 if ((iand(rdr(iwrd+i*10+10),flagd).eq.0.or.(flagw.and.iand(rdr(iwrd+i*10+10),32768).ne.0))!weird bit
     & .and.spot(i)
     & .and.rdr(iwrd+i*10+1).ne.invalid !geolocation
     & .and.rdr(iwrd+i*10+4).ne.-1      !range 
     & ) then !no nonsense
	  dlat(i) = dscale*rdr(iwrd+i*10+2)
	  dlon(i) = dscale*rdr(iwrd+i*10+1)
	  if(geographic .and. dlon(i).lt.0.) dlon(i)=dlon(i)+360.d0
          rnge(i) = rscale*unsigned(rdr(iwrd+i*10+4))
	  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
c decide output data
	  if(rflag) then
	   if(rdr(iwrd+i*10+4) .eq. -1) goto 10
	   data(i) = rnge(i)
	  else if(flagp) then
	   if(rdr(iwrd+i*10+5) .eq. -1) goto 10 !pulsewidths do not pass go, do not collect $200!
	   data(i) = pscale*rdr(iwrd+i*10+5)
	  else if(bidf) then !bidirectional reflectivity
c albedo
c trans(i) contains aperture, transmission efficiency, etc
c scale factor accounts for km^2, mJ to fJ
	   enrg(i) = escale*rdr(iwrd+i*10+6)
	   data(i)=rnge(i)**2 *escale*(enrg(i)-dark(i))/(xenrg*trans(i))
	  else !default is height
	   data(i) = rscale*unsigned(rdr(iwrd+i*10+3)) - ref !height
	  endif
	  irec1(i)=irec1(i)+1
	  write(10+i,rec=irec1(i))dlon(i),dlat(i),data(i)
	  if(extend.and.dlon(i).ge.(dlonmin+360.d0)) then
	   irec1(i)=irec1(i)+1
	   write(10+i,rec=irec1(i)) dlon(i)-360.d0,dlat(i),data(i)
	  endif
	  if(extend.and.dlon(i).le.(dlonmax-360.d0)) then
	   irec1(i)=irec1(i)+1
	   write(10+i,rec=irec1(i)) dlon(i)+360.d0,dlat(i),data(i)
	  endif
 10	  continue
	 endif! spot flags
c	 nois(i) =        rdr(iwrd+i*10+7)
c	 thrs(i) = tscale*rdr(iwrd+i*10+8)
c	 gain(i) = gscale*rdr(iwrd+i*10+9)
c	 gflg(i) = iand(1,rdr(iwrd+i*10+10))
	 endif ! lonlat
	enddo
	endif ! limits
        goto 100 !big loop
999     continue
        if(spot(0).and.orec.gt.0) then
	 close(10,status='keep')
	else
	 close(10,status='delete')
	endif
	do i=1,5
	 if(irec1(i).gt.0) then
	  disp='keep  '
	  write(*,*)'spot ',i,' rec=',irec1(i)
	 else
	  disp='delete'
	 endif
	 close(unit=i+10,status=disp)
	enddo
	goto 99999
9999	write(0,*)'file not opened',filename
99999   end