c program pedr2tab.f for (small-endian) PC using MSFortran
c ----------------------------------------------------------------------
c output ASCII table of MOLA shot values from PEDR binary files
c The PDS label is checked for SFDU and software keyword.
c ----------------------------------------------------------------------
c input preferences from PEDR2TAB.PRM
c ----------------------------------------------------------------------
c Version 2.7, applies parallax and crossover correction
c Refer to PEDR Software Interface Specification, Version 2.7,
c http://ltpwww.gsfc.nasa.gov/tharsis/pedrsis-2-7.html
c
c ----------------------------------------------------------------------
c written in Fortran-77 for IEEE-std, big-endian machines, with some
c language extensions and system functions that may cause errors:
c
c variable names may be longer than, but are unique to six letters
c
c The $ format descriptor extension is used to output to the
c format buffer without end-of-line terminators.
c
c The byte data type extension is used to decode signed 1 byte values,
c sometimes referred to as integer*1
c
c Direct access I/O is used to read binary unformatted data.
c For DEC compilers, redefine the "recl=776" keyword to "recl=194"
c
c ----------------------------------------------------------------------
USE MSFLIB
implicit real*8 (a-h,o-z)
parameter (pi=3.141592653589793d0)
parameter (d2r=pi/180.d0,r2d=180.d0/pi)
parameter (irecl=776)!bytes for all but DEC compilers
c parameter (irecl=194)!words for DEC VAX,ALPHA
integer*4 iswap4, a
integer*2 iswap2, j,iargloc
logical lhdr, lall,lgrd,lpr(0:7),obf,usgs, first,lxovr,list
character*64 filelist, file2, filestring, filetab, filename
character*16 ltmn,ltmx,lnmn,lnmx
character*8 ccsd
real*4 yth(4)
integer bkgrd,mask
integer*2 iwant,loctime,iphase,incidence,iemission,aflag
integer*2 mgmver,mgm2
c The PEDR data frame:
integer*4 pedr(194)
integer*2 ipdr(388)
byte bpdr(776)
equivalence (pedr,ipdr,bpdr)
c this equivalence for string search in header:
character*776 cpdr
equivalence (cpdr,bpdr)
c some values are accessed conveniently through equivalences:
equivalence (loctime, ipdr(271))
equivalence (iphase, ipdr(272))
equivalence (incidence, ipdr(273))
equivalence (iemission, ipdr(274))
data filelist /'FILELIST'/
data file2 /'MOLA.TAB'/
data dlatdr /0./, dlondr/0./, xcorr /0./,dxcor/0./!versions <3
data dltmn/-90.d0/,dltmx/90.d0/,dlnmn/-360.d0/,dlnmx/360.d0/
c
c data mask/#3FFF/ !packet counter bits 13-0
data mask/z'3FFF'/ !packet counter bits 13-0
data yth /2.29,1.32,0.763,0.440/ !threshold voltage gain factors
data lall /.FALSE./ ! all shot returns
data lgrd /.TRUE./ ! ground returns
data iwant /1/ ! shot_classification_flag =0 => clouds or noise,
c =1 => ground return
data lhdr /.TRUE./,obf /.FALSE./! two header lines, one big file
data lpr /1*.TRUE.,7*.FALSE./, first /.TRUE./! selected values
data usgs /.FALSE./ ! areographic coordinates?
data lxovr /.TRUE./ ! crossover corrections if true
data list /.TRUE./ ! read from filelist
data mgm2 /2/ !first version with parallax
c PEDR 4-byte integer values:
c 1: frame mid-point whole seconds past J2000
c 2: frame mid-point fractional microseconds past J2000
c 3: rev number
c 4: sc areocentric lat deg*1000000
c 5: sc lon deg*1000000
c 6: mgs radial distance
c 7: mid-point range (cm)
c 8: bad shot bits (1=bad)
c 9-12 quality bit flags
c 13-32 shot planetary radii (cm)
c 33: midpoint planetary radius (cm)
c 34: right ascension
c 35: declination
c 36: twist
c --> 37-46 rcv pulse energy
c --> 47-56 reflectivity-transmissivity product
c --> 57-61 channel number bytes
c --> 62-71 pulsewidth at threshold crossing
c --> 72-81 pulsewidth at threshold crossing
c 82-83: parallax dlat/dr and dlon/dr
c 84: drift function derived from crossovers
c 85: frame midlat deg*1000000
c 86: frame midlon deg*1000000
c --> 87-96 laser transmit energy
c --> 97-106 shot classification codes
c --> 107-114 background noise counts
c 115 range window delay, cm
c 116 range window width, cm
c --->117-134 other engineering data
c --->135-137: flags and angles
c 138: opacity
c 139-140: frame midpoint time, IEEE-double
c 140-150: pulse width-area counts
c 151: delta_sc_lat
c 152: delta_sc_lon
c 153: delta_sc_rad cm
c 154 areoid radius
c 155 off-nadir pointing angle, deg*1000000
c 161 delta areoid
c 162 MOLA clock frequency
c 163-182 MOLA_range (cm)
c 193 delta-latitude
c 194 delta-longitude
c ipdr 2-byte integer values:
c 246 frame index, 1-7
c 248 packet number, bits 13-0
c 254 fine time code
c 269 orbit code
c 270 attitude code
c 271 local time
c 272 solar phase
c 273 solar incidence
c 274 laser emission
c 365-384 range correction (cm)
c bpdr 1-byte values
c 225-244: trigger channel
c 561-580: received energy counts
c 581-600: pulse width counts
c ----------------------------------------------------------------------
iswap2(j)=ior(ishft(iand(j,#FF),8),iand(ishft(j,-8),#FF))
c
iswap4(a)=ior(ior( ior(ishft(a,24),iand(ishft(a,8),#00FF0000)),
1 iand(ishft(a,-8),#0000FF00)),
2 iand(ishft(a,-24),#000000FF))
c ----------------------------------------------------------------------
c execution begins with flattening default value (Viking MDIMS):
aob2 = (3393.4d0/3375.73d0)**2
c input format preferences
open(unit=1,file='PEDR2TAB.PRM',status='old',err=1)
read(1,*,end=1,err=1) lhdr
do i=0,7
read(1,*,end=1,err=1) lpr(i)
enddo
usgs = lpr(2) !Areographic needed
read(1,*,end=1,err=1) lall
read(1,*,end=1,err=1) lgrd
read(1,*,end=1,err=1) lxovr
read(1,*,end=1,err=1) obf , file2
c write(*,*)'File2=', file2
read(1,*,end=1,err=1) dlnmn
read(1,*,end=1,err=1) dlnmx
read(1,*,end=1,err=1) dltmn
read(1,*,end=1,err=1) dltmx
read(1,*,end=1,err=1) flatn
aob2 = 1./(1. - 1./flatn)**2
c read(1,*,end=1,err=1) usgs
goto 2
1 continue
write(*,*)' PEDR2TAB.PRM not found or wrong format'
2 continue
close(1)
c ----------------------------------------------------------------------
write(*,*) ' Output parameters (PEDR2TAB.PRM): '
write(*,*) ' -------------------------------- '
write(*,*) ' Header lines: ',lhdr
write(*,*) ' shot location, topo, flags: ',lpr(0)
write(*,*) ' MGS location: ',lpr(1)
write(*,*) ' angle, ET, areodetic_lat, areoid: ',lpr(2)
write(*,*) ' shot #, packet #, rev #, GMM #: ',lpr(3)
write(*,*) ' solar time, phase, incidence: ',lpr(4)
write(*,*) ' range walk & pulse statistics: ',lpr(5)
write(*,*) ' background, threshold, raw pulse: ',lpr(6)
write(*,*) ' range window, range delay: ',lpr(7)
write(*,*) ' selected (F) or all (T) shots: ',lall
write(*,*) ' noise/clouds (F), ground shots (T):',lgrd
if(lgrd) then
iwant =1
else
iwant =0
endif
write(*,*) ' apply crossover corrections: ',lxovr
write(*,*) ' template filetype or one big file :',obf
id=index(file2,".")
write(*,10)' lon. (positive East):',dlnmn,' to',dlnmx
write(*,10)' lat. (areocentric): ',dltmn,' to',dltmx
10 format(2(a,f9.2))
if (usgs) write(*,10) ' Flattening (areographic): ',flatn
c ----------------------------------------------------------------------
if (iargc().gt.0) then
iargloc=1
call getarg(iargloc,filelist)
else
write(*,*)
write(*,*)' PEDR or List of Binary PEDR Files:'
read(*,'(a)') filelist
endif
open(unit=4,file=filelist,status='old',action='read',err=9999)
if(index(filelist,".b").gt.0 .or.index(filelist,".B").gt.0)then
list = .false.!assume PEDR, check first 8 chars of SFDU
read(4,'(a)',end=999) ccsd
if(ccsd.eq."CCSD3ZF0") then
filestring=filelist
close(4)
else
write(*,*)'Invalid PEDR header! ', ccsd
goto 9999
endif
else !leave 4 open for loop
write(*,*)'List of PEDR files'
list = .true.
endif
c One Big File?
if(obf) open(unit=2,file=file2,status='unknown',
& iostat=iopen,err=9999)
ict=0
c ----------------------------------------------------------------------
c Loop over input files
c ----------------------------------------------------------------------
4 continue
if (list) read(4,'(a)',end=99) filestring
idot=index(filestring,".b")
if (idot.eq.0) idot=index(filestring,".B")
c write(*,*) idot, filestring(1:idot)
if (idot.eq.0) then
write(*,*)'Filename must have valid extension'
write(*,*)'skipping "',filestring,'"'
close (1)
if (list) goto 4
endif
write(*,'(a,a)')' Opening:',filestring(1:idot+1)
open(unit=1,file=filestring,status='old',iostat=iopen,
& access='direct',action='read',recl=irecl,err=99)
inquire(unit=1, name=filename)
idot=index(filename,".")
c test version
read(1,rec=1,err=9) pedr !label
isoft=index(cpdr,"SOFTWARE_NAME =")
if (isoft.gt.0) then
read(unit=cpdr(isoft+25:isoft+28),fmt=*) version
write(*,'(a,f5.2)')' software is',version
if (version.lt. 7.125) then
write(*,*) cpdr(isoft:isoft+29),' lacks xovers'
lxovr = .false.
endif
if (version.le.7.0d0) mgm2 = 3
else
c die?
write(*,*) 'SOFTWARE_NAME not found'
if (list) goto 4
endif
c ----------------------------------------------------------------------
if(.not. obf) then
filetab = filename(1:idot)//file2(id+1:id+3)
open(unit=2,file=filetab,
& status='unknown',
& iostat=iopen,err=999)
write(*,*)'writing: ', filetab
first = .true.
endif
if (lhdr .and. first) then
c headings
if(lpr(0)) write(2,8000)
if(lpr(1)) write(2,8001)
if(lpr(2)) write(2,8002)
if(lpr(3)) write(2,8003)
if(lpr(4)) write(2,8004)
if(lpr(5)) write(2,8005)
if(lpr(6)) write(2,8006)
if(lpr(7)) write(2,8007)
write(2,8008)
if(lpr(0)) write(2,8010)
if(lpr(1)) write(2,8011)
if(lpr(2)) write(2,8012)
if(lpr(3)) write(2,8013)
if(lpr(4)) write(2,8014)
if(lpr(5)) write(2,8015)
if(lpr(6)) write(2,8016)
if(lpr(7)) write(2,8017)
write(2,8008)
first = .false.
endif
read(1,rec=11,err=9) pedr !first frame
dptime=iswap4(pedr(1)) + 1.d-06*iswap4(pedr(2))
etime0=(dptime-1.) !time of first frame
c BAG to set slope of correction function at second frame
icorr=iswap4(pedr(84))
mgmver=iswap2(ipdr(269))
if (mgmver.ge.mgm2) then
c crossover corrector, once per frame
lastxc = icorr
read(1,rec=12,err=9) pedr
icorr=iswap4(pedr(84))
lastxc = lastxc - (icorr-lastxc)
endif
c
ifound=0
do i=1,99999
read(1,rec=i+10,err=9) pedr
dptime=iswap4(pedr(1)) + 1.d-06*iswap4(pedr(2))
irev = iswap4(pedr(3))
iseq=iand(iswap2(ipdr(248)),mask)
dlat= 1.d-6*iswap4(pedr(85))
dlon= 1.d-6*iswap4(pedr(86))
c Check frame midpoint for bounds
if (dlat.ge.dltmn.and.dlat.le.dltmx
& .and.dlon.ge.dlnmn.and.dlon.le.dlnmx) then
ifound=ifound+1
mgmver=iswap2(ipdr(269))
c parallax and drift corrections
if (mgmver.ge.mgm2) then
dlatdr =1.d-9*iswap4(pedr(82))
dlondr =1.d-9*iswap4(pedr(83))
c crossover corrector, once per frame
if (lxovr) then
icorr=iswap4(pedr(84))
xcorr = 1.d-2*icorr
dxcor = 1.d-2*(icorr-lastxc)
lastxc = icorr
endif
endif
c mola clock frequency, scaled
f = 1.d-9*iswap4(pedr(162))
flat= 1.d-6*iswap4(pedr(4))
flon= 1.d-6*iswap4(pedr(5))
fscrad=1.d-2*iswap4(pedr(6))
dmidrange=1.d-2*iswap4(pedr(7))
fplrad=1.d-2*iswap4(pedr(33))
rdelay=1.d-2*iswap4(pedr(115))
rwnd=1.d-2*iswap4(pedr(116))
dslatdt=1.d-6*iswap4(pedr(151))
dslondt=1.d-6*iswap4(pedr(152))
draddt=1.d-2*iswap4(pedr(153))
fref =1.d-2*iswap4(pedr(154))
drefdt=1.d-2*iswap4(pedr(161))
dlatdt=1.d-6*iswap4(pedr(193))
dlondt=1.d-6*iswap4(pedr(194))
aflag=iswap2(ipdr(270))
c scale angles
hourlocal = 12.d0 + iswap2(loctime)*0.0012d0/pi
solarp = iswap2(iphase)*0.0180d0/pi
solari = iswap2(incidence)*0.0180d0/pi
emission= iswap2(iemission)*0.0180d0/pi
offndr = 1.d-6*iswap4(pedr(155))
do k=1,20
molarg=iswap4(pedr(k+162))
if (molarg.gt.0 ! non-zero range returned
& .and. (lall .or.iswap2(ipdr(k+192)).eq. iwant)) then !output shot
ict=ict+1
c for each half-frame
i2 = (k-1)/10
ichan=bpdr(k+224)
c shot planetary radius (raw)
plrad = 0.01d0*iswap4(pedr(k+12))
x=(k-10.5d0)/20.d0
c MGS location
flatk = flat + x*dslatdt
flonk = flon + x*dslondt
if(flonk.lt.0.) flonk=flonk+360.d0
if(flonk.ge.360.d0) flonk=flonk-360.d0
c shot location corrected for parallax
dlatk = dlat + x*dlatdt+dlatdr*(plrad-fplrad)
dlonk = dlon + x*dlondt+dlondr*(plrad-fplrad)
if(dlonk.lt.0.) dlonk=dlonk+360.d0
if(dlonk.ge.360.d0) dlonk=dlonk-360.d0
c altimetry
plrad = plrad-xcorr-x*dxcor! after crossover corrections
range = 0.01d0*molarg
areoid = fref + x*drefdt
topo = plrad - areoid
dmgsrad=fscrad+x*draddt
etimex = dptime +(0.01d0/f)*(k-10.5d0)
if(usgs) phi= r2d*atan(aob2*tan(d2r*dlatk))
Rc = 0.01d0*iswap2(ipdr(k+364))
iframe = iswap2(ipdr(246))
ishot=k+20*(iframe-1) ! 1-140 shot count
c "corrected, scaled received_pulse_energy"
ErfJoule=iswap2(ipdr(k+72))
ireft = iswap2(ipdr(k+92))
if(ireft.lt.0) ireft=ireft+65536
Reft=0.001*ireft ! Percent Reflectivity*Transmissivity
c raw counts stored as byte values
ipact=bpdr(k+560)
if(ipact.lt.0) ipact=ipact+256
ipwct=bpdr(k+580)
Pwt=0.1*iswap2(ipdr(k+122))
Sigopt=0.1*iswap2(ipdr(k+142))
Elaser=0.01*iswap2(ipdr(k+172))
thrsh=iswap2(ipdr(232+4*i2+ichan))*yth(ichan)
c plog2 values are converted already in pproc
bkgrd=iswap4(pedr(106+4*i2+ichan))
c long lat topo range plan.radius chan att_flag
if(lpr(0))
& write(2,9000)dlonk,dlatk,topo,range,plrad,ichan,aflag
if(lpr(1))write(2,9001) flonk, flatk, dmgsrad
if(lpr(2))write(2,9002) offndr,etimex,phi,areoid
if(lpr(3))write(2,9003) ishot,iseq,irev,mgmver
if(lpr(4))write(2,9004) hourlocal,solarp,solari
if(lpr(5))write(2,9005)
& emission,Rc,Pwt,Sigopt,Elaser,ErfJoule,Reft
if(lpr(6))write(2,9006) bkgrd,thrsh,ipwct,ipact
if(lpr(7))write(2,9007) rwnd,rdelay
write(2,8008)
endif
enddo !end of shot loop
endif !end of bounds check
c
enddo !end of i/o loop
9 continue
close(unit=1)
write(*,'(a,a)')' Closing: ', filetab
if(.not.obf) close(unit=2)
write(*,*)' Frames read: ', ifound
write(*,*)' Total shots output: ', ict
if (list) goto 4
99 continue
close(unit=2)
999 write(*,*) ' Done!'
close(4)
c ----------------------------------------------------------------------
c format statements
c ----------------------------------------------------------------------
9000 format(f9.5,f10.5,f11.2,f10.2,f12.2,2i2,$)
9001 format(2f10.5, f12.2,$)
9002 format(f7.3,f16.5,f10.5,f11.2,$)
9003 format(i4,i6,i5,i3,$)
9004 format(f7.3,2f7.2,$)
9005 format(f7.2,f7.2,2f8.1,f8.2,f7.0,f7.3,$)
9006 format(i7, f7.1,i3,i4,$)
9007 format(f7.0,f8.0,$)
8000 format(
& 'long_East lat_North topography MOLArange planet_rad c A',$)
8010 format(
& '---.----- --.----- ------.-- ------.-- --------.-- - -',$)
8001 format(
& ' SC_long SC_lat SC_radius',$)
8011 format(
& ' ---.----- ---.----- --------.--',$)
8002 format(
& ' offndr EphemerisTime areod_lat areoid_rad',$)
c8002 format(
c & ' offndr EphemerisTime areoid_rad',$)
8012 format(
& ' --.--- --------------- --.----- ------.--',$)
c8012 format(
c & ' --.--- --------------- ------.--',$)
8003 format(
& ' shot pkt orbit gm',$)
8013 format(
& ' --- ---- ---- --',$)
8004 format(
& ' hrlcl s_phas s_inc ',$)
8014 format(
& ' --.--- ---.-- ---.--',$)
8005 format(
& ' emissn Rcorr PWT Sigopt Elaser PulseE Ref*T%',$)
8015 format(
& ' ---.-- ----- ---.- ------ ----- ------ --.---',$)
8006 format(
& ' bkgrd th_mV Wct Ect',$)
8016 format(
& ' ----- ------ -- ---',$)
8007 format(
& ' R_wnd R_dly ',$)
8017 format(
& ' ------ -------',$)
8008 format(a)
c ----------------------------------------------------------------------
c normal return or error handling
c ----------------------------------------------------------------------
goto 99999
9999 continue
write(*,*)'Error ',iopen,' opening files'
99999 end