-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathRinexSNRv2.f
237 lines (225 loc) · 8.46 KB
/
RinexSNRv2.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
program RinexSNRv2
implicit none
include 'local.inc'
integer stderr
parameter (stderr=6)
character*80 inline
character*80 rawfilename, outfilename, broadfile
character*4 station
character*2 prn_pickc
character*1 satID(maxsat)
integer nobs, itime(5), prn(maxsat), numsat, flag, sec,
+ msec, lli(maxob,maxsat), iprn, l2c_pts,
. ios, itrack, L2Ctracking, i, iobs(maxsat), jtime(5),
. gpsweek, the_hour, prn_pick, current_hour,npts,
. blockIIRM(maxsat)/maxsat*0/, fileIN, fileOUT,iymd(3),
. rinex_year, rinex_doy,inside_year,inside_doy,nav_doy,
. nav_year, ts, k1, k2
real*8 obs(maxob,maxsat), tod,North(3),
. East(3), Up(3), azimuth, elev, staXYZ(3), tod_save,
. s1,s2,xrec, yrec, zrec, tc, l1,l2, Lat,Long,Ht,
. bele(maxeph, maxsat, 28),grange,s5,
. nxrec, nyrec, nzrec, tjul, rho,tjul_jan1
logical eof, bad_point
logical help, fsite, debugging
c Kristine M. Larson
c version 2.0 new version - uses subroutines, fixed bug in az-el
c verison 2.1 fixes bug in selection 77 (reading the LLI value)
c version 2.2 fixed bug in reading the comment lines, 15sep07
c version 2.3 added S5
c version 2.5 distribution on gitHub
c version 2.6 added moving site velocities, station marker name
c version 2.7 fixed bug for reading files with < 5 observations
c version 2.8 changed minimum observables from 10 to 15
c version 2.9 19mar01 - allow up to 20 observables, changed code
c to read the header and the data block
c version 2.10 check dates on input files to make sure they are
c in agreement
c
c common block for the GPS broadcast ephemeris
include 'new_orbit.inc'
c define the blockIIRM satellites
call pickup_blockIIRM(blockIIRM)
c set some defaults - used for input and output file IDs
debugging = .false.
l2c_pts = 0
fileIN = 12
fileOUT = 55
tod_save = 0.0
npts = 0
bad_point = .false.
help = .false.
c read input files
call getarg (1,rawfilename)
call getarg (2,outfilename)
call getarg (3,broadfile)
call getarg (4,prn_pickc)
write(stderr, *) '----------------------------------'
write(stderr, *) 'Input RINEX: ', rawfilename
write(stderr, *) 'Broadcast RINEX: ', broadfile
write(stderr, *) 'Outputfile: ', outfilename
if (rawfilename(1:2) .eq. ' ') help = .true.
c
if (help) then
write(stderr,*)
+ 'gpsSNR.e inputRinex outputSNR navigation-file out-choice '
write(stderr,*) ' '
write(stderr,*) 'out-choice can be a specific PRN number or'
write(stderr,*) '50 - all satellites, all data < 10 degrees'
write(stderr,*) '66 - all satellites, all data < 30 degrees'
write(stderr,*) '77 - L2C satellites, all data > 5 degrees'
write(stderr,*) '88 - all satellites, all data > 5 degrees'
write(stderr,*) '99 - all satellites, 5 < data < 30 degrees'
call exit
endif
c figure out which option is being requested
READ (prn_pickc, '(I2)') prn_pick
write(stderr, *) 'Selection ', prn_pick
c Check to see if broadcast file exists
open(22,file=broadfile, status='old',iostat=ios)
if (ios.ne.0) then
print*, 'Problem opening navigation file:',broadfile
print*, 'Exiting'
close(22)
call exit(0)
endif
close(22)
c read the header of the RINEX file, returning station coordinates
c 19mar01 - change to 20 observables
call read_header_20obs(fileIN,rawfilename, xrec,yrec,zrec,
. iobs,nobs,iymd,station)
call name2ydoy(rawfilename,rinex_year, rinex_doy)
call name2ydoy(broadfile,nav_year, nav_doy)
write(stderr,*)'The year in your filename : ', rinex_year
write(stderr,*)'The day of year in your filename: ', rinex_doy
call check_dates(iymd,rinex_year,rinex_doy,
. inside_year,inside_doy,nav_year,nav_doy)
c KL added 18jan15
call moving_sites(station, iymd(1), iymd(2), iymd(3),
. nxrec,nyrec,nzrec,fsite)
if (fsite) then
print*, 'Using variable model station coordinates (meters)'
xrec = nxrec
yrec = nyrec
zrec = nzrec
print*, xrec, yrec, zrec
endif
if ((xrec+yrec+zrec) .eq.0) then
print*, 'No apriori coordinates - exiting'
call exit
endif
if (iobs(6) .eq. 0) then
print*, 'no L1 SNR data - exiting '
call exit
endif
c use the proper limit instead of hardwiring 20
if (nobs .gt. maxobs) then
print*, 'This code only works for <='
print*, maxob, ' observations types'
call exit
endif
c read the broadcast ephemeris information
call read_broadcast4 (broadfile, bele,iymd)
call rearrange_bele(bele)
call envTrans(xrec,yrec,zrec,staXYZ,Lat,Long,Ht,North,East,Up)
c open output file
open(fileOUT,file=outfilename, status='unknown')
eof = .false.
current_hour = 0
print*, 'S1 location:', iobs(6)
print*, 'S2 location:', iobs(7)
print*, 'S5 location:', iobs(8)
c start reading the observation records
do while (.not.eof)
inline = ' '
read(fileIN,'(A80)', iostat=ios) inline
if (ios.ne.0) goto 99
read(inline(1:32),'(5I3,X,I2,X,I3,4X,2I3)')
+ (itime(i), i=1,5), sec, msec, flag, numsat
if (itime(4) .ne.current_hour) then
current_hour = itime(4)
endif
c seconds in the day
tod = itime(4)*3600.0 + 60.0*itime(5) + sec
if (tod.lt.tod_save) then
c print*,'Time is going backwards or standing still'
c print*,'Ignoring this record'
bad_point = .true.
tod_save = tod
else
bad_point = .false.
tod_save = tod
endif
if (debugging) then
print*, 'reading block '
print*, inline(1:60)
endif
c 19mar01 - expanding number of observables allowed
call read_block_gps(fileIN, flag,inline,numsat,nobs,satID,
. prn,obs,lli)
c if flag has value 4, that means there were comment
c lines, and those were skipped
if (flag .ne. 4) then
call convert_time(itime,sec, msec, gpsweek, tc)
do itrack = 1, numsat
c only uses GPS satellites currently, skips over other
c constellations
if ((satID(itrack).eq.'G').or.(satID(itrack).eq.' ')) then
the_hour = 1+nint(1.d0*itime(4) + itime(5)/60)
if (the_hour > 24) the_hour = 24
iprn = prn(itrack)
c get azimuth and elevation angle
call get_azel(tc, iprn, staXYZ,East,North,Up,
. Lat,Long,Ht,azimuth,elev,the_hour,tod,grange)
c you can modify this to also pick up pseudoranges
c l1 = obs(iobs(1),itrack)
c l2 = obs(iobs(2),itrack)
c these can be modfiied to store phas data. currently 0
l1 = 0.0
l2 = 0.0
s1 = obs(iobs(6),itrack)
s2 = obs(iobs(7),itrack)
c check for S5
if (iobs(8).ne.0) then
s5 = obs(iobs(8),itrack)
else
s5=0
endif
c check to see if L2C is available in the file
if (iobs(2) .ne. 0) then
L2Ctracking = lli(iobs(2),itrack)
else
L2Ctracking = 0
endif
if (debugging) then
print*, s1, s2, s5
endif
if (s1.eq.0.d0 .and. s2.eq.0.d0) then
c no data, SNR so do not print it
else
call write_to_file(fileOUT, prn_pick, iprn,
. elev,azimuth, tod, s1,s2, bad_point,l1,l2,
. grange,blockIIRM, L2Ctracking,s5)
if ((prn_pick .eq. 77).and.(s2.gt.0.d0)) then
l2c_pts = l2c_pts + 1
endif
npts = npts + 1
endif
else
c print*, 'skipping non-GPS satellite'
endif
enddo
endif
enddo
99 continue
if (npts.eq.0) then
write(stderr,*)
. 'You have been misled. There are no S1/S2 data in this file'
endif
c close input and output files
if (prn_pick .eq. 77) then
print*, 'L2C points ', l2c_pts
endif
close(fileIN)
close(fileOUT)
end