program g12rin *** convert garmin gps 12xl log file into rinex implicit double precision(a-h,o-z) character*6 verdat logical unstuf dimension nmess(0:255) byte b1,b2,b3,b4,b5,b6,buff(255+4),bx byte DLE,ETX parameter(DLE = z'10') parameter(ETX = z'03') logical lhdr save /misc/ common/misc/lrnx,igwk,lhdr *** version date verdat='020203' lin=1 lrnx=2 open(lin, file='garbin.bin',form='binary',status='old') open(lrnx,file='g12xl.rnx',form='formatted',status='unknown') 1 write(*,'(a$)') ' GPS week = ' read(*,*,err=1) igwk if(igwk.le.0) go to 1 call rinini call rinhdr(verdat) do 2 i=1,255 2 nmess(i)=0 *** *** loop over messages until eof *** 100 read(lin,end=777) b1,b2,b3,b4 20 read(lin,end=777) b5 if(b1.ne.DLE.or.b2.ne.ETX.or.b3.ne.DLE) then b1=b2 b2=b3 b3=b4 b4=b5 go to 20 endif messid=mybyte(b4) nmess(messid)=nmess(messid)+1 meslen=mybyte(b5) if(.not.unstuf(meslen,lin,bx)) then if(bx.eq.ETX) then b1=meslen b2=bx read(lin,end=777) b3,b4 go to 20 endif endif *** data packet do 30 i=1,meslen read(lin,end=777) buff(i) if(.not.unstuf(buff(i),lin,bx)) then !*** short record if(bx.eq.ETX) then b1=buff(i) b2=bx read(lin,end=777) b3,b4 go to 20 endif endif 30 continue read(lin,end=777) b6 !*** checksum if(.not.unstuf(b6,lin,bx)) then !*** short record if(bx.eq.ETX) then b1=b6 b2=bx read(lin,end=777) b3,b4 go to 20 endif endif if(messid.eq.z'38') call do38(buff) go to 100 *** *** end of file encountered *** 777 call dumprb do 3 i=0,255 if(nmess(i).gt.0) * write(*,'(a,z3.2,a,i8)') 'id=',i,' count=',nmess(i) 3 continue end logical function unstuf(b,lin,b2) *** unstuff dle stuffing byte DLE,b,b2 parameter(DLE = z'10') if(b.eq.DLE) then read(lin,end=777) b2 if(b2.ne.DLE) then unstuf=.false. return endif endif unstuf=.true. return 777 stop 98765 end subroutine do38(buff) *** process message 38 (measurement data block) implicit double precision(a-h,o-z) parameter(maxsat=12,iepflg=0) parameter(maxprn=32) parameter(sol= 299792458.d0) parameter(F1 = 1575.42d6) dimension iprns(maxsat),phases(maxsat),ll1s(maxsat), * isnrs(maxsat),prs(maxsat) logical dec5,lstlok,oppsgn integer getQcode byte buff(*) integer*4 phscnt byte tbyte byte b1,b2,b3 integer*2 deltaf integer*4 intphs real*8 pr integer*4 c511500 integer*2 isigq real*8 tow byte bp real*8 ba2dbl integer*4 ba2int integer*2 ba2in2 phscnt = ba2int(buff( 1)) tbyte = buff( 5) *** skip bytes 6-8 deltaf = ba2in2(buff( 9)) intphs = ba2int(buff(11)) pr = ba2dbl(buff(15)) c511500 = ba2int(buff(23)) isigq = ba2in2(buff(27)) tow = ba2dbl(buff(29)) bp = buff(37) *** 5 second interval if(mod(idnint(tow),5).ne.0) return iprn = mybyte(bp)+1 icfrac = iand(phscnt,2#00000000000000000000011111111111) phase = dble(intphs) + icfrac/2048.d0 isnr = getQcode(isigq) *** From RINEX document: *** Time(corr) = Time(r) - dT(r) *** PR(corr) = PR(r) - dT(r)*c *** phase(corr) = phase(r) - dT(r)*freq dt = tow - idnint(tow) tow = tow - dt pr = pr - dt*sol phase = phase - dt*F1 *** check tracking and status if(tbyte.eq.0) return if(iprn .le.0) return if(iprn .gt.maxprn) return if(tow .le.1.d-9) return if(tow .gt.1604801.d0) return if(pr .le.1400000.d0) return if(pr .gt.30000000.d0) return if(isigq.le.0) return call putrin(tow,iprn,pr,phase,isnr) return *** end of file encountered 777 stop end *********************************************************************** double precision function ba2dbl(b) *** convert byte array to double byte b(*),ba(8) double precision db equivalence(ba,db) ba(1)=b(1) ba(2)=b(2) ba(3)=b(3) ba(4)=b(4) ba(5)=b(5) ba(6)=b(6) ba(7)=b(7) ba(8)=b(8) ba2dbl=db return end integer function ba2int(b) *** convert byte array to integer byte b(*),ba(4) integer in equivalence(ba,in) ba(1)=b(1) ba(2)=b(2) ba(3)=b(3) ba(4)=b(4) ba2int=in return end integer*2 function ba2in2(b) *** convert byte array to short byte b(*),ba(2) integer*2 in2 equivalence(ba,in2) ba(1)=b(1) ba(2)=b(2) ba2in2=in2 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 getQcode(sgn) *** convert internal signal strength to standardized value *** note: static method -- use for all instances integer*2 sgn if (sgn.ge.9000) then getQcode=9 return elseif(sgn.ge.6000) then getQcode=8 return elseif(sgn.ge.4000) then getQcode=7 return elseif(sgn.ge.2500) then getQcode=6 return endif getQcode=5 return end subroutine rinhdr(verdat) *** rinex header implicit double precision(a-h,o-z) character*6 verdat integer*2 iyr2,imo2,idy2,ihr2,imn2,isc2,ish2 logical lhdr save /misc/ common/misc/lrn,igwk,lhdr 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)') * 'g12rin ',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)') * '0002 ', * 'MARKER NUMBER' write(lrn,'(a,a)') * 'your name here ', * 'OBSERVER / AGENCY' write(lrn,'(a,a)') * 'nnnnnnnn GARMIN 12XL 4.57 ', * 'REC # / TYPE / VERS' write(lrn,'(a,a)') * 'GARMIN 12XL Internal Antenna ', * '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)') * ' 2 0 ', * 'WAVELENGTH FACT L1/2' write(lrn,'(a,a)') * ' 2 C1 L1 ', * '# / TYPES OF OBSERV' write(lrn,'(a,a)') * ' 5.0000 ', * 'INTERVAL' write(lrn,'(a,a)') * 'Forced Modulo Decimation to 5 seconds ', * 'COMMENT' write(lrn,'(a,a)') * ' SNR is mapped to RINEX snr flag value [1-9] ', * 'COMMENT' *** enable header writing second part lhdr=.true. return end subroutine rinhd2(towx) *** rinex header implicit double precision(a-h,o-z) logical lhdr save /misc/ common/misc/lrnx,igwkx,lhdr igwk=igwkx tow =towx call sowciv(igwk,tow,iyr,imo,idy,ihr,imn,sec) write(lrnx,'(5i6.2,f13.7,a,a)') * iyr,imo,idy,ihr,imn,sec,' GPS ', * 'TIME OF FIRST OBS' write(lrnx,'(a,a)') * ' ', * 'END OF HEADER' *** disable any more header writing lhdr=.false. return end *********************************************************************** *** rinex data buffer *** subroutine rinini *** initialize rinex buffer implicit double precision(a-h,o-z) parameter(mst=12,maxprn=32) save /rinbuf/ common/rinbuf/prs(mst),phss(mst),iprns(mst),isnrs(mst),tow,nsat nsat=0 tow=-1.d0 return end subroutine putrin(towx,iprn,pr,phs,isnr) *** put rinex data into buffer (support for rinex epoch header) implicit double precision(a-h,o-z) logical dupprn parameter(mst=12) parameter(ttol=0.001) !*** tolerance (sec.) for tow control save /rinbuf/ common/rinbuf/prs(mst),phss(mst),iprns(mst),isnrs(mst),tow,nsat *** dump if new tow -- else just load if(dabs(tow-towx).ge.ttol) then if(tow.ge.0.d0) then call dumprb else call rinhd2(towx) endif tow=towx nsat=1 prs(nsat) =pr phss(nsat) =phs iprns(nsat)=iprn isnrs(nsat)=isnr *** increment and load buffer else if(nsat.lt.mst) then if(.not.dupprn(iprn)) then nsat=nsat+1 prs(nsat) =pr phss(nsat) =phs iprns(nsat)=iprn isnrs(nsat)=isnr endif else write(*,*) 'rinex buffer overflow' stop 23843 endif endif return end logical function dupprn(iprn) *** test for duplicate prn in one epoch implicit double precision(a-h,o-z) parameter(mst=12,maxprn=32) save /rinbuf/ common/rinbuf/prs(mst),phss(mst),iprns(mst),isnrs(mst),tow,nsat do 1 i=1,nsat if(iprn.eq.iprns(i)) then dupprn=.true. return endif 1 continue dupprn=.false. return end subroutine dumprb *** dump the rinex buffer implicit double precision(a-h,o-z) parameter(mst=12) save /rinbuf/ common/rinbuf/prs(mst),phss(mst),iprns(mst),isnrs(mst),towx,nsat save /misc/ common/misc/lrnx,igwkx *** epoch header first igwk=igwkx tow =towx call sowciv(igwk,tow,iyr,imo,idy,ihr,imn,sec) if(iyr.gt.1999) then !*** y2k iyr=iyr-2000 else iyr=iyr-1900 endif write(lrnx,1) iyr,imo,idy,ihr,imn,sec,nsat,(iprns(i),i=1,nsat) 1 format(5i3,f11.7,' 0',i3,12('G',i2)) *** first pseudorange, then phase do 10 i=1,nsat 10 write(lrnx,2) prs(i),isnrs(i),phss(i),isnrs(i) 2 format(5(f14.3,' ',i1)) 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