program ant2rin *** convert u-blox antaris ubx (binary) log file into rinex *** done ------------------------------------ *** new snr map *** drop bad phase option *** always drop waas *** 5 sec decimation (manage phase flags) *** Message QI option *** bug fix -- insure no epoch header if no sats implicit double precision(a-h,o-z) character*6 verdat logical lhdr dimension nmess(0:255,0:255) integer*2 mleng byte b1,b2,b3,b4,buff(1255+4) byte bclass,bmesid *** version date verdat='060509' lin=1 lrnx=2 open(lin, file='ant2rin.ubx',form='binary',status='old') open(lrnx,file='ant2rin.obs',form='formatted',status='unknown') call doargs call inislip call rinhdr(lrnx,verdat,lhdr) do i=0,255 do j=0,255 nmess(i,j)=0 !*** initialize message counters enddo enddo call slipset !*** set all cycle slip flags false *** *** loop over messages until eof *** 100 read(lin,end=777) b1 200 read(lin,end=777) b2 if(mybyte(b1).ne.181.or.b2.ne.98) then b1=b2 go to 200 endif *** synch characters detected -- read prefix read(lin,end=777) bclass,bmesid,mleng mclass=mybyte (bclass) messid=mybyte (bmesid) meslen=myshort(mleng) *** process payload (or skip message contents) if(mclass.eq.2.and.messid.eq.16) then call doraw(lin,lrnx,lhdr) else read(lin,end=777) (buff(i),i=1,meslen) endif if(mclass.lt.0.or.mclass.gt.255) stop 39493 if(messid.lt.0.or.messid.gt.255) stop 39494 nmess(mclass,messid)=nmess(mclass,messid)+1 !*** count message *** skip the checksums read(lin,end=777) b3,b4 !*** 2 checksum bytes go to 100 *** end of file encountered 777 continue *** message count do i=0,255 do j=0,255 if(nmess(i,j).gt.0) * write(*,'(a,2i4,a,i9)') 'cls,id=',i,j,' count=',nmess(i,j) enddo enddo *** write slip counts call getnslip(nmqi4,nmqi5,nmqi6) write(*,'(a,i6)') '# MQI=4 drops: ',nmqi4 write(*,'(a,i6)') '# MQI=5 drops: ',nmqi5 write(*,'(a,i6)') '# MQI=6 drops: ',nmqi6 end *********************************************************************** subroutine doraw(lin,lrnx,lhdr) *** process RXM-RAW (measurement measurement data) payload implicit double precision(a-h,o-z) parameter(mxsat=16,iepflg=0) dimension iprns(mxsat),phases(mxsat),prs(mxsat),dops(mxsat), * mqis(mxsat),cnos(mxsat),llis(mxsat),isnrs(mxsat),idrp(mxsat) logical lhdr,ldec,lbit0,lbit1 logical lqual,lfull,ldrop common/args/lqual,lfull,ldrop,igran save /args/ parameter(mxprn=32) logical lslip(0:mxprn),loppf(0:mxprn) common/llock/lslip,loppf save /llock/ *** note: SVN is really PRN code integer*4 itow integer*2 igweek integer*1 mesqi,isnr byte bnsv,bdumm1,bprn,blli real*8 cpm,prm real*4 dom *** read header read(lin) itow,igweek,bnsv,bdumm1 nsv=mybyte(bnsv) if(nsv.lt.0.or.nsv.gt.mxsat) return !*** abort this record nsv2=nsv !*** number of output satellites *** gps time in seconds of week igwk=igweek !*** signed I2 into signed I4 gpst=itow/1000.d0 !*** milliseconds to seconds call sowciv(igwk,gpst,iyr,imo,idy,ihr,imn,sec) *** decimation flag (if ldec=.false. -- do not output epoch) if(igran.le.0) then ldec=.true. else if(mod(idnint(gpst*10.d0),igran*10).eq.0) then !*** igran (s.) ldec=.true. else ldec=.false. endif endif *** loop over satellite blocks -- read and parse do i = 1,nsv read(lin) cpm,prm,dom,bprn,mesqi,isnr,blli phases(i)=cpm !*** L1 cyc. prs (i)=prm !*** m dops (i)=dom !*** Hz iprns (i)=mybyte(bprn) mqis (i)=mesqi !*** ubx-rxm-raw message quality cnos (i)=dble(isnr) !*** dbHz llis (i)=mybyte(blli) !*** rinex loss-of-lock indicator isnrs (i)=mapsnr(isnr) idrp (i)=1 !*** initialize flag (do not drop) if(iprns(i).gt.mxprn) then !*** drop waas (set iprn=0) iprns(i)=0 idrp (i)=0 !*** set flag to drop nsv2=nsv2-1 !*** decrement # kept satellites endif ***** sticky bit for cycle slip history if( .not.lslip(iprns(i)) .and. lbit0(llis(i)) ) * lslip(iprns(i))=.true. ***** use most recent opposite factor bit (do not keep history) loppf(iprns(i))=lbit1(llis(i)) enddo *** scan data -- drop sats with bad phase (set idrp=0) *** only pass the best phases if(ldrop) then do i = 1,nsv if(dabs(phases(i)).le.0.0005d0.or.mqis(i).lt.7) then if (mqis(i).eq.6) then call incr6 elseif(mqis(i).eq.5) then call incr5 else call incr4 endif if(idrp(i).ne.0) then idrp(i)=0 !*** set flag to drop nsv2=nsv2-1 !*** decrement # kept satellites endif endif enddo endif *** output epoch header and measurement list *** lli (bit 0 or 1 set is a cycle slip, bit 2 is AS on) *** snr is mapped [1-9] *** don't output epoch if decimation flag not set *** don't output sat if idrp was set to 0 (drop flag) *** skip header if no sats to be written if(ldec) then if(nsv2.gt.0) then if(lhdr) call rinhd2(lrnx,iyr,imo,idy,ihr,imn,sec,lhdr) write(lrnx,'(5i3,f11.7,i3,i3$)') * mod(iyr,100),imo,idy,ihr,imn,sec,iepflg,nsv2 do i = 1,nsv if(idrp(i).ne.0) then if(iprns(i).le.0.or.iprns(i).gt.mxprn) stop 98484 write(lrnx,'(''G'',i2$)') iprns(i) endif enddo write(lrnx,'(a1)') endif if(lqual) then do i = 1,nsv if(idrp(i).ne.0) then llis(i)=lliset(lslip(iprns(i)),loppf(iprns(i))) write(lrnx,'(f14.3,2i1,f14.3,2i1,f14.3,2i1,f14.3,2i1)') * phases(i),llis(i),isnrs(i), * prs (i),llis(i),isnrs(i), * dble(mqis (i)),llis(i),isnrs(i), * cnos (i),llis(i),isnrs(i) endif enddo elseif(lfull) then do i = 1,nsv if(idrp(i).ne.0) then llis(i)=lliset(lslip(iprns(i)),loppf(iprns(i))) write(lrnx,'(f14.3,2i1,f14.3,2i1,f14.3,2i1,f14.3,2i1)') * phases(i),llis(i),isnrs(i), * prs (i),llis(i),isnrs(i), * dops (i),llis(i),isnrs(i), * cnos (i),llis(i),isnrs(i) endif enddo else do i = 1,nsv if(idrp(i).ne.0) then llis(i)=lliset(lslip(iprns(i)),loppf(iprns(i))) write(lrnx,'(f14.3,2i1,f14.3,1x,i1)') * phases(i),llis(i),isnrs(i),prs(i),isnrs(i) endif enddo endif call slipset !*** set all cycle slip flags false endif !*** ldec return end *********************************************************************** *** Loss of lock indicator (LLI). Range: 0-7 *** 0 or blank: OK or not known *** Bit 0 set : Lost lock between previous and *** current observation: cycle slip *** possible *** Bit 1 set : Opposite wavelength factor to the *** one defined for the satellite by a *** previous WAVELENGTH FACT L1/2 line. *** Valid for the current epoch only. *** Bit 2 set : Observation under Antispoofing *** (may suffer from increased noise) *** *** Bits 0 and 1 for phase only. *********************************************************************** subroutine rinhdr(lrn,verdat,lhdr) *** rinex header implicit double precision(a-h,o-z) character*6 verdat logical lhdr integer*2 iyr2,imo2,idy2,ihr2,imn2,isc2,ish2 logical lqual,lfull,ldrop common/args/lqual,lfull,ldrop,igran save /args/ call getdat(iyr2,imo2,idy2) call gettim(ihr2,imn2,isc2,ish2) write(lrn,'(a,a)') * ' 2.10 OBSERVATION DATA G (GPS) ', * 'RINEX VERSION / TYPE' write(lrn,'(a,a,a,i4,2i2.2,i3.2,a,i2.2,a,i2.2 ,a,a)') * 'ant2rin ',verdat,' your name here ', * iyr2,imo2,idy2,ihr2,':',imn2,':',isc2,'ET ', * 'PGM / RUN BY / DATE' write(lrn,'(a,a)') * 'XXXX ', * 'MARKER NAME' write(lrn,'(a,a)') * '0001 ', * 'MARKER NUMBER' write(lrn,'(a,a)') * 'your name here ', * 'OBSERVER / AGENCY' write(lrn,'(a,a)') * 'nnnnn u-blox LEA-4T Firmware V5.00 ', * 'REC # / TYPE / VERS' write(lrn,'(a,a)') * '600577 u-blox ANN-MS-0-005 ', * 'ANT # / TYPE' write(lrn,'(a,a)') * ' 1095790.2358 -4831326.5883 4003934.2615 ', * 'APPROX POSITION XYZ' write(lrn,'(a,a)') * ' 0.0000 0.0000 0.0000 ', * 'ANTENNA: DELTA H/E/N' write(lrn,'(a,a)') * ' 1 0 ', * 'WAVELENGTH FACT L1/2' if(lqual) then write(lrn,'(a,a)') * ' 4 L1 C1 D1 S1 ', * '# / TYPES OF OBSERV' write(lrn,'(a,a)') * 'Doppler field really contains MessQI (float format) ', * 'COMMENT' elseif(lfull) then write(lrn,'(a,a)') * ' 4 L1 C1 D1 S1 ', * '# / TYPES OF OBSERV' else write(lrn,'(a,a)') * ' 2 L1 C1 ', * '# / TYPES OF OBSERV' endif ***debug -- making assumptions on data rates if(igran.eq.5) then write(lrn,'(a,a)') * ' 5.0000 ', * 'INTERVAL' write(lrn,'(a,a)') * 'Forced Modulo Decimation to 5 seconds ', * 'COMMENT' else write(lrn,'(a,a)') * ' 1.0000 ', * 'INTERVAL' endif write(lrn,'(a,a)') * ' SNR is mapped to RINEX snr flag value [1-9] ', * 'COMMENT' write(lrn,'(a,a)') * ' L1: 10 dBHz = 1; 30 dBHz = 5; 50 dBHz = 9 ', * 'COMMENT' *** enable header writing second part lhdr=.true. return end subroutine rinhd2(lrn,iyr,imo,idy,ihr,imn,sec,lhdr) *** rinex header implicit double precision(a-h,o-z) logical lhdr write(lrn,'(5i6.2,f13.7,a,a)') * iyr,imo,idy,ihr,imn,sec,' GPS ', * 'TIME OF FIRST OBS' write(lrn,'(a,a)') * ' ', * 'END OF HEADER' *** disable any more header writing lhdr=.false. return end integer function mapsnr(isnr) *** map snr to rinex interval *** 10 dBHz = 1; 30 dBHz = 5; 50 dBHz = 9 implicit double precision(a-h,o-z) mapsnr = idnint((isnr-10)/5.d0)+1 if(mapsnr.lt.1) mapsnr=1 if(mapsnr.gt.9) mapsnr=9 return end *********************************************************************** subroutine inislip *** initialize slip counts implicit integer(a-z) common/pslps/nmqi4,nmqi5,nmqi6 save /pslps/ nmqi4=0 nmqi5=0 nmqi6=0 return end subroutine getnslip(mqi4,mqi5,mqi6) *** return the slip counts implicit integer(a-z) common/pslps/nmqi4,nmqi5,nmqi6 save /pslps/ mqi4=nmqi4 mqi5=nmqi5 mqi6=nmqi6 return end subroutine incr4 *** increment the MessQI 4 count implicit integer(a-z) common/pslps/nmqi4,nmqi5,nmqi6 save /pslps/ nmqi4=nmqi4+1 return end subroutine incr5 *** increment the MessQI 5 count implicit integer(a-z) common/pslps/nmqi4,nmqi5,nmqi6 save /pslps/ nmqi5=nmqi5+1 return end subroutine incr6 *** increment the MessQI 6 count implicit integer(a-z) common/pslps/nmqi4,nmqi5,nmqi6 save /pslps/ nmqi6=nmqi6+1 return end *********************************************************************** subroutine doargs *** process command line arguments *** caution: system-dependent fortran 77 (MS-fort) implicit double precision(a-h,o-z) character*255 buff integer*2 ipos,status n=nargs() !*** including the program name do i=2,n !*** skip the program name ipos=i-1 call getarg(ipos,buff,status) istat=status if(istat.gt.0) call doarg(buff,istat) enddo return end subroutine iniargs *** initialize default argument values implicit double precision(a-h,o-z) logical lqual,lfull,ldrop common/args/lqual,lfull,ldrop,igran save /args/ lqual=.false. !*** do not dump qual.i. & s.stength lfull=.false. !*** do not dump doppler & s.stength ldrop=.false. !*** do not drop phase igran=0 !*** do not decimate data return end subroutine doarg(buff,istat) *** initialize default argument values implicit double precision(a-h,o-z) character*(*) buff logical lqual,lfull,ldrop common/args/lqual,lfull,ldrop,igran save /args/ *** -q, quality indicator dump if(buff(1:2).eq.'-q') then write(*,*) 'q -- quality indicator dump' lqual=.true. lfull=.false. *** -f, full dump (disable if quality indicator dump) elseif(buff(1:istat).eq.'-f'.and. .not.lqual) then write(*,*) 'f -- full dump of data' lfull=.true. *** -d, drop bad phase elseif(buff(1:2).eq.'-d') then write(*,*) 'd -- drop bad phase' ldrop=.true. *** -r, decimation to 5 seconds elseif(buff(1:2).eq.'-r') then write(*,*) 'r -- decimate to 5 seconds' igran=5 endif return end *********************************************************************** subroutine slipset *** set all slip flags to false implicit double precision(a-h,o-z) parameter(mxprn=32) logical lslip(0:mxprn),loppf(0:mxprn) common/llock/lslip,loppf save /llock/ do i=0,mxprn lslip(i)=.false. loppf(i)=.false. enddo return end logical function lbit0(lli) *** is bit 0 set on integer*4 loss of lock indicator? *** caution: system-dependent fortran 77 (MS-fort) logical btest lbit0 = btest(lli,0) return end logical function lbit1(lli) *** is bit 1 set on integer*4 loss of lock indicator? *** caution: system-dependent fortran 77 (MS-fort) logical btest lbit1 = btest(lli,1) return end integer function lliset(lslp,lopp) *** construct a loss of lock indicator *** assume no A/S for CA code-derived carrier phase ***--------------------------------------------------- *** Loss of lock indicator (LLI). Range: 0-7 *** 0 or blank: OK or not known *** Bit 0 set : Lost lock between previous and *** current observation: cycle slip *** possible *** Bit 1 set : Opposite wavelength factor to the *** one defined for the satellite by a *** previous WAVELENGTH FACT L1/2 line. *** Valid for the current epoch only. *** Bit 2 set : Observation under Antispoofing *** (may suffer from increased noise) *** *** Bits 0 and 1 for phase only. ***--------------------------------------------------- implicit integer(a-z) logical lslp,lopp lliset=0 if(lslp) lliset=lliset+1 if(lopp) lliset=lliset+2 return end *********************************************************************** integer function mybyte(b) *** convert byte to integer byte b if(b.lt.0) then mybyte=b+256 else mybyte=b endif return end integer function myshort(iu2) *** convert short unsigned integer to integer integer*2 iu2 if(iu2.lt.0) then myshort=iu2+65536 else myshort=iu2 endif return end subroutine sowciv(igwk,sow,iyr,imo,idy,ihr,imn,sec) *** convert gps week/seconds of week to civil implicit double precision(a-h,o-z) parameter(mjd6jan80=44244) *** fmjdg computation preserves digits mjd = igwk*7 + sow/86400.d0 + mjd6jan80 fmjd = dmod(sow,86400.d0)/86400.d0 call mjdciv(mjd,fmjd,iyr,imo,idy,ihr,imn,sec) return end subroutine mjdciv(mjd,fmjd,iyr,imo,idy,ihr,imn,sec) *** convert modified julian date to civil date *** imo in range 1-12, idy in range 1-31 *** only valid in range mar-1900 thru feb-2100 *** ref: hofmann-wellenhof, 2nd ed., pg 34-35 *** operation confirmed for leap years (incl. year 2000) implicit double precision(a-h,o-z) rjd=mjd+fmjd+2400000.5d0 ia=(rjd+0.5d0) ib=ia+1537 ic=(ib-122.1d0)/365.25d0 id=365.25d0*ic ie=(ib-id)/30.6001d0 *** the fractional part of a julian day is fractional mjd + 0.5 *** therefore, fractional part of julian day + 0.5 is fractional mjd it1=ie*30.6001d0 idy=ib-id-it1+fmjd it2=ie/14.d0 imo=ie-1-12*it2 it3=(7+imo)/10.d0 iyr=ic-4715-it3 tmp=fmjd*24.d0 ihr=tmp tmp=(tmp-ihr)*60.d0 imn=tmp sec=(tmp-imn)*60.d0 return end