program spcsnr * This program will fit SNR values and estimate corrections * to phase data based on the oscillations of the SNR. include 'spcsnr.h' * ierr - IOSTAT error * rcpar - Reads the runstring * len_run - Length of the runstring read * i,j - Loop counters * iter - Number of iterattions integer*4 ierr, rcpar, len_run, i,j, iter * more_data - Logical to indicate more sequences of data available * flat - Logical to indicate that the spectrum is flat * split_seq - Logical to indicate that sequence has been divided in * half logical more_data, flat, split_seq character*32 runstring **** Get the name of the input snr file (generated by svsnr) len_run = rcpar(1, infile) if( len_run.le.0 ) then call proper_runstring('spcsnr.hlp','spcsnr',-1) end if open(100,file=infile, status='old', iostat=ierr) call report_error('IOSTAT',ierr,'open',infile,1,'spcsnr') **** Tell user what is happen write(*,120) infile(1:len_run) 120 format('* SPCSNR: Spectral fitting SNR data in ',a) * Get minium elevation len_run = rcpar(2, runstring ) if( len_run.gt.0 ) then read(runstring,*, iostat = ierr) min_elev call report_error('IOSTAT',ierr,'Min ELEV decod',runstring, . 1,'SPCSNR') else * Set default 10 degrees if none passed or bad * value passed. min_elev = 10.d0 end if * Get the minum SNR values len_run = rcpar(3, runstring ) if( len_run.gt.0 ) then read(runstring,*, iostat = ierr) min_snr call report_error('IOSTAT',ierr,'Min SNR decod',runstring, . 1,'SPCSNR') else * Set default value of 2 for Min SNR * value passed. min_snr = 2.d0 end if **** Decode the output options len_run = rcpar(4, runstring ) * Replace the separators blanks so that the call sub_char(runstring,':',' ' ) * set default output to 40 == DPHS and RMSE only out_opts = 40 call decode_option(runstring, out_cmds, 32, out_opts, 40) **** See if input c-file is given to be updated in_cf = ' ' len_run = rcpar(5, in_cf ) **** Check for update cffile name out_cf = ' ' **** Now read all of the input file call read_snr(100) write(*,220) num_epochs, num_sat 220 format('* ',i5,' epochs of data with ',i3,' satellites found') * Now fit the elevation angle dependent model to the SNR data num_poly_L1 = max_poly num_poly_L2 = max_poly call fit_elev * Clear the phase adjusments do i = 1, num_epochs do j = 1, num_sat phs_l1a(i,j) = 0.0 phs_L2a(i,j) = 0.0 end do end do **** Now loop over each satellite to get the corrections do i = 1, num_sat * Start by pulling off the segments of data for this * satellite more_data = .true. split_seq = .false. start_seq = 0 do while ( more_data ) * Get the next segment of data call ext_data( i, more_data, split_seq, 'L1') * Now Process L1 first: flat = .false. num_cmp = 0 iter = 0 do while ( .not.flat .and. more_data .and. . num_seq.gt.10 .and. iter.lt.50 ) iter = iter + 1 * Copy the current instrument gains into the data * arrays for FFT. call copy_fft( gain_L1q, seq, num_seq, max_seq) call fourg( seq, num_seq, 1, work, 6) * now find the maxium spectral peak call find_max_spec( seq, num_seq, mf, flat, out_opts ) * If the sequence is not flat yet, then do the * fit for the components. This routine also removes * the spectral components from gain sequence so that * next clean component can be removed. if( .not.flat ) then call est_snr_spec( gain_L1q ) end if * Compute the phase corrections call sum_dphs( gain_L1s, phs_L1a, i) end do * Update the sequence positions start_seq = end_seq end do * Now do the L2 Gain results. more_data = .true. split_seq = .false. start_seq = 0 do while ( more_data ) * Get the next segment of data call ext_data( i, more_data, split_seq, 'L2' ) flat = .false. num_cmp = 0 iter = 0 do while ( .not.flat .and. more_data .and. . num_seq.gt.10 .and. iter.lt.50 ) iter = iter + 1 * Copy the current instrument gains into the data * arrays for FFT. call copy_fft( gain_L2q, seq, num_seq, max_seq) call fourg( seq, num_seq, 1, work, 6) * now find the maxium spectral peak call find_max_spec( seq, num_seq, mf, flat, out_opts ) * If the sequence is not flat yet, then do the * fit for the components. This routine also removes * the spectral components from gain sequence so that * next clean component can be removed. if( .not.flat ) then call est_snr_spec( gain_L2q ) end if * Compute the phase corrections call sum_dphs( gain_L2s, phs_L2a, i) end do * Update the sequence positions start_seq = end_seq end do end do **** Now write out the final phase correction values call out_dphs end CTITLE COPYR subroutine copyr * This routine copies the difference bwteen the smooth and * observed SNR values back into the observed arrays so that * they can be smoothed less tightly. include 'spcsnr.h' * LOCAL VARIABLES * i,j,k - Loop counters integer*4 i,j * kbit - Checks the bits in a word logical kbit ***** Loop doing the smoothing do i = 1, num_epochs do j = 1, num_sat if( kbit(flags(i,j),1) ) then gain_L1o(i,j) = 1.d0 + gain_L1o(i,j) - gain_L1s(i,j) gain_L2o(i,j) = 1.d0 + gain_L2o(i,j) - gain_L2s(i,j) end if end do end do ***** Thats all return end CTITLE OUT_PHS subroutine out_dphs * This routine writes out the adjustments to the phase values * (L1 and L2 are found from the SNR, LC and LG are computed * from these. include 'spcsnr.h' include 'const_param.h' * LOCAL VARIABLES * i,j,k - Loop counters integer*4 i,j, date(5), code_type * dl1, dl2, dlc, dlg - Adjustments to L1, L2, LC and LG real*8 dl1, dl2, dlc, dlg, sectag, sod * Statistics for elevation statistics real*8 sum_el(18), var_el(18), mean, rms, . sum_al, var_al integer*4 bin, num_el(18), ne, num_al * kbit - Checks the bits in a word logical kbit ***** Clear statistics do i = 1,18 num_el(i) = 0 sum_el(i) = 0 var_el(i) = 0 end do num_al = 0 sum_al = 0 var_al = 0 ***** Loop doing the smoothing if( kbit(out_opts,4) ) then write(*,110) min_elev, min_snr 110 format('* Phase corrections computed with ',F5.2, . ' deg minimum elevation, and ',F5.2, . ' Min SNR DPHS',/, . '* Date h m s sod PRN Azimuth', . ' Elev L1 L2 LC LG C', . ' DPHS'/, . '* (deg) ', . ' (deg) (cyc) (cyc) (cyc) (cyc) ', . ' DPHS') endif do i = 1, num_epochs call mjd_to_ymdhms( epoch(i)+1.d-8, date, sectag) sod = date(4)*3600.d0 + date(5)*60.d0 + sectag do j = 1, num_sat if( kbit(flags(i,j),1) ) then dL1 = phs_L1a(i,j) dL2 = phs_L2a(i,j) dLC = (dL1 - (fL2/fL1)*dL2)/(1-(fL2/fL1)**2) dLG = -(dL1-(fL1/fL2)*dL2)/(1.d0-(fL1/fL2)**2) code_type = 1 if( kbit(flags(i,j),3) ) code_type = 0 if( kbit(out_opts,4) ) . write(*,120) date, sectag, sod, prn_list(j), . az(i,j)*180/pi, el(i,j)*180/pi, . dL1, dL2, dLC, dLG, code_type 120 format(1x,i4,4i3,f6.2,1x,f8.2,' PRN ',i2.2,2f9.4, . 4f9.5,1x,i2, ' DPHS ') * Accumulate statistics bin = (el(i,j)*180/pi)/5+1 if( bin.gt.18 ) bin = 18 if( bin.lt.1 ) bin = 1 if( dLC.ne.0.d0 ) then num_el(bin) = num_el(bin) + 1 sum_el(bin) = sum_el(bin) + dLC var_el(bin) = var_el(bin) + dLC**2 num_al = num_al + 1 sum_al = sum_al + dLC var_al = var_al + dLC**2 end if end if end do end do ***** Write out the Elevation statistics if( kbit(out_opts,6) ) then write(*,220) 220 format(': Elevation angle statistics for LC',/ . ': Elev Num Mean (cyc) RMS (cyc)') do i = 1,18 ne = num_el(i) if( ne.gt. 0 ) then mean = sum_el(i)/ne rms = sqrt(var_el(i)/ne) write(*,230) (i-0.5)*5.0, ne, mean, rms 230 format(':',f6.2,1x,i6,2f10.3) end if end do if( num_al.gt.0 ) then mean = sum_al/num_al rms = sqrt(var_al/num_al) write(*,240) num_al, mean, rms 240 format(': ALL ',1x,i6,2f10.3) end if end if ***** Write out the smoothed gains if( kbit(out_opts,5) ) then write(*,310) 310 format('* Smoothed Normalized SNR from FFT fit SMTH',/, . '* Date h m s sod PRN ', . ' Azimuth Elev L1 NSNR L2 NSNR C SMTH') do i = 1, num_epochs call mjd_to_ymdhms( epoch(i)+1.d-8, date, sectag) sod = date(4)*3600.d0 + date(5)*60.d0 + sectag do j = 1, num_sat if( kbit(flags(i,j),1) ) then code_type = 1 if( kbit(flags(i,j),3) ) code_type = 0 write(*,320) date, sectag, sod, prn_list(j), . az(i,j)*180/pi, el(i,j)*180/pi, . gain_L1s(i,j), gain_L2s(i,j), code_type 320 format(1x,i4,4i3,f6.2,1x,f8.2,' PRN ',i2.2,2f9.4, . 2f10.5,1x,i2, ' SMTH ') end if end do end do end if ! Outputing the Smoothed Normalized SNR. ***** Thats all return end CTITLE READ_SNR subroutine read_snr( unit ) * Routine to read the SNR file generated by svsnr. include 'spcsnr.h' include 'const_param.h' * PASSED Variables * unit - Unit number integer*4 unit * LOCAL Variables * ierr, jerr - IOSTAT errors. * i,j - Loop counter * trimlen - Length of string * chan - Channel number for the prn read * prn - PRN number * date(5) - Date of measurement * code_type - Indicates if X-correlation (1) or code * correlation (0) integer*4 ierr, jerr, i, trimlen, chan, prn, date(5), . code_type * sectag - Seconds tage * sod - Seconds of day (not used) * jdread - Julian date computed from read date * lastjd - Last jd read * azread, elread - Azimith and elevation angles read (degrees) * snr1, snr2 - Signal-to-noise ratios read from file real*8 sectag, sod, jdread, lastjd, azread, elread, . snr1, snr2 * line - Line read from file character*256 line **** Start reading the file num_epochs = 0 num_sat = 0 num_L2_code(1) = 0 num_L2_code(2) = 0 ierr = 0 do while ( ierr.eq.0 ) read(unit,'(a)', iostat=ierr ) line if( ierr.eq.0 .and. line(1:1).eq.' ' .and. * ! Process the line . trimlen(line).gt.0 ) then read(line,120,iostat=jerr) date, sectag, sod, prn, . azread, elread, snr1, snr2, code_type 120 format(i5,4i3,F6.1,f10.1,1x,'PRN ',i2.2,1x, . 2f10.3,1x,2f6.1,3x,i2) * Remove negative SNR if( snr1.le.0 ) snr1 = 0 if( snr2.le.0 ) snr2 = 0 call ymdhms_to_mjd( date, sectag, jdread) **** See if the epoch has changed if( abs(jdread-lastjd).gt.1.d0/86400.d0 ) then num_epochs = num_epochs + 1 lastjd = jdread epoch(num_epochs) = jdread do i = 1, max_sat flags(num_epochs,i) = 0 end do end if **** Now see what channel the Prn is in call get_chan( prn, prn_list, num_sat, chan) **** OK save the values in the appropriate locations az(num_epochs,chan) = azread*pi/180.d0 el(num_epochs,chan) = elread*pi/180.d0 * Check the SNR value to see if above limit if( snr1.lt.min_snr ) snr1 = 0.d0 snr_L1o(num_epochs,chan) = snr1 * Remove the L1 contribution from the snr if this * is not code tracking (ie. X-correlation) if( snr1.eq.0 .or. snr2.lt.min_snr) snr2 = 0.d0 if( code_type.ne.0 ) then snr_L2o(num_epochs,chan) = 100.d0*snr2/snr1 num_L2_code(2) = num_L2_code(2) + 1 else snr_L2o(num_epochs,chan) = snr2 num_L2_code(1) = num_L2_code(1) + 1 end if * Set flags bit to show that we have an SNR measurement * at this time call sbit(flags(num_epochs,chan),1,1) * Set L1 bit to show that it is code measurement call sbit(flags(num_epochs,chan),2,1) * If L1 SNR greter than 0 then mark OK if( snr1.gt.0 .and. elread.gt.min_elev) then call sbit(flags(num_epochs,chan),4,1) end if * If L2 SNR greter than 0 then mark OK if( snr2.gt.0 .and. elread.gt.min_elev) then call sbit(flags(num_epochs,chan),5,1) end if * Check to see if code measurement at L2, if it is * set bit to show this. (If not a code measurement then * L1 SNR gain will be removed before fiting the L2 gain curve). if( code_type.eq.0 ) then call sbit(flags(num_epochs,chan),3,1) end if end if end do **** Thats all return end CTITLE GET_CHAN subroutine get_chan( prn, prn_list, num_sat, chan ) * This routine will return the chan number for a given PRN * or add the PRN to the list * prn - PRN number to be matched * prn_list(*) - List of PRN's already found * num_sat - Number of PRN's already found * chan - The returned channel number integer*4 prn, prn_list(*), num_sat, chan * LOCAL Variables * i - Loop counter integer*4 i **** See if we ready have i = 1 do while ( prn_list(i).ne.prn .and. i.le.num_sat ) i = i + 1 end do if ( prn_list(i).ne. prn ) then num_sat = num_sat + 1 prn_list(num_sat) = prn chan = num_sat else chan = i end if **** Thats all return end CTITLE FIT_ELEV subroutine fit_elev * This routine will fit a polynomial in sin(elev) to the L1 and * L2 SNR values * The basic L1 model is * SNR_L1 = A sin(el) + B sin(el)**2 + .. * For L2 (X-correlation) we first remove the L1 SNR model (scale to 1 * at zenith) and then fit a second set of coefficients to sin(el). For * the L2 code tracking channels an over-all scale factor is estimated. include 'spcsnr.h' include 'const_param.h' * LOCAL Variables * i,j,k - Loop counters * ipivot(max_poly) - Pivot elements for inversion * date(5) - Date of measurement * code_type - Code type (set to 0 for X-correlation, 1 for code * - correlation) integer*4 i,j,k, ipivot(max_poly), date(5), code_type * kbit - Checks bit status logical kbit * se - Sin(elevation) * apart(max_poly) - Paritial derivatives for estimates * scale(max_poly) - Scale factors for inversion * est_snr - Estimated SNR from polynomial fit * sectag - Seconds tag of date * sod - Seconds of day. * elv - Elevatiob * PG_L1, PG_L2 - Computed gains from L1 and L2 from Polynomial fit. real*8 se, apart(max_poly), scale(max_poly), est_snr, . sectag, sod, elv, PG_L1, PG_L2 * cov_col(max_poly,1) - Covariance matrix colum for forcing * sol_col(1) - Solution variable * dchi - Change in chi**2 * avat(1,1) - AvAT matrix * equ_gn(max_poly,1) - Gain matrix real*8 cov_col(max_poly,1), sol_col(1), dchi, avat(1,1), . equ_gn(max_poly,1) * nc(1) - Parameter numbers to be forced. integer*4 nc(1) **** Start with the L1 SNR fit. call clear_norm(num_poly_L1) * Now loop over all of the data. do i = 1, num_epochs do j = 1, num_sat * ! SNR good if( kbit(flags(i,j),1) .and. kbit(flags(i,j),4) ) then se = sin(el(i,j)) do k = 1, num_poly_L1 apart(k) = se**(k-1) end do **** Now increment the normal equations call inc_norm(snr_l1o(i,j), apart, num_poly_L1) end if end do end do **** Now solve the system of equations call invert_vis( norm_eq, bvec, scale, ipivot,num_poly_L1, . max_poly, 1 ) **** OK, Save the gain curve values zen_L1 = 0 do i = 1, num_poly_L1 gain_poly_L1(i) = bvec(i) zen_L1 = zen_L1 + bvec(i) end do write(*,120) zen_L1, (gain_poly_L1(i),i-1, i = 1, num_poly_L1) 120 format('* L1 Zenith SNR: ',F8.2,/, . '* L1 Gain curve: ',10(F8.2,'*sin(el)**',i1,1x)) write(*,125) num_L2_code 125 format('* There are ',i8,' L2 code measurements, ',i8, . ' cross correlation measurements') **** See if we should force the constant term to zero because it is less than * zero. if( gain_poly_L1(1).lt.0 ) then nc(1) = 1 call force_parms( norm_eq, bvec, max_poly, cov_col, sol_col, . nc, 1, 0.d0, dchi, avat, equ_gn, ipivot, . scale ) zen_L1 = 0 do i = 1, num_poly_L1 gain_poly_L1(i) = bvec(i) zen_L1 = zen_L1 + bvec(i) end do write(*,120) zen_L1, (gain_poly_L1(i),i-1, i = 1, num_poly_L1) end if ***** Now compute the observed L1 gains, and set up L2 by removing L1 gain do i = 1, num_epochs do j = 1, num_sat * ! SNR good if( kbit(flags(i,j),1) ) then se = sin(el(i,j)) est_snr = 0 * Exclude the zero-order term from gain. Directly * subtract this from the SNR. This seems to be needed * for Astech-Z12 data. do k = 1, num_poly_L1 est_snr = est_snr + gain_poly_L1(k)*se**(k-1) end do gain_L1o(i,j) = snr_L1o(i,j)/est_snr if( gain_L1o(i,j).lt.0 ) then gain_L1o(i,j) = 0.d0 call sbit(flags(i,j),4,0 ) end if * * Save the L2 gain gain_L2o(i,j) = snr_L2o(i,j) end if end do end do **** Now fit the L2 elevation dependence call clear_norm(num_poly_L2) * Now loop over all of the data. do i = 1, num_epochs do j = 1, num_sat * ! SNR good and if( kbit(flags(i,j),5) ) then se = sin(el(i,j)) do k = 1, num_poly_L2 apart(k) = se**(k-1) end do * Accumulate solution depending on whether * there are more or less X-correlation measurements. if( num_L2_code(1).lt.num_L2_code(2) .and. . .not. kbit(flags(i,j),3) ) then * ! X-correlation **** Now increment the normal equations call inc_norm(gain_L2o(i,j), apart, num_poly_L2) else if ( num_L2_code(1).ge.num_L2_code(2) .and. . kbit(flags(i,j),3) ) then * ! code measurement **** Now increment the normal equations call inc_norm(gain_L2o(i,j), apart, num_poly_L2) end if end if end do end do **** Now solve the system of equations call invert_vis( norm_eq, bvec, scale, ipivot,num_poly_L2, . max_poly, 1 ) **** OK, Save the gain curve values zen_L2 = 0 do i = 1, num_poly_L2 gain_poly_L2(i) = bvec(i) zen_L2 = zen_L2 + bvec(i) end do write(*,220) zen_L2,(gain_poly_L2(i),i-1, i = 1, num_poly_L2) 220 format('* L2 Zenith SNR: ',F8.2,/, . '* L2 Gain curve: ',10(F8.2,'*sin(el)**',i1,1x)) **** See if we should force the constant term to zero because it is less than * zero. if( gain_poly_L2(1).lt.0 ) then nc(1) = 1 call force_parms( norm_eq, bvec, max_poly, cov_col, sol_col, . nc, 1, 0.d0, dchi, avat, equ_gn, ipivot, . scale ) zen_L2 = 0 do i = 1, num_poly_L2 gain_poly_L2(i) = bvec(i) zen_L2 = zen_L2 + bvec(i) end do write(*,220) zen_L2, (gain_poly_L2(i),i-1, i = 1, num_poly_L2) end if ***** Now compute the observed L2 gains, and set up L2 by removing L1 gain norm_eq(1,1) = 0 bvec(1) = 0 do i = 1, num_epochs do j = 1, num_sat * ! SNR good if( kbit(flags(i,j),1) ) then se = sin(el(i,j)) est_snr = 0 * Remove zero-order term do k = 1, num_poly_L2 est_snr = est_snr + gain_poly_L2(k)*se**(k-1) end do ***** Update gain if X-correlated and fit to X-correlation * of if code and fit to code. gain_L2o(i,j) = gain_L2o(i,j)/ est_snr if ( gain_l2o(i,j).lt.0 ) then gain_l2o(i,j) = 0 call sbit(flags(i,j),5,0 ) end if if( num_L2_code(1).lt. num_L2_code(2) .and. . kbit(flags(i,j),3) .and. kbit(flags(i,j),5) ) then norm_eq(1,1) = norm_eq(1,1) + 1 bvec(1) = bvec(1) + gain_L2o(i,j) else if ( num_L2_code(1).ge. num_L2_code(2) .and. . .not.kbit(flags(i,j),3) .and. . kbit(flags(i,j),5) ) then norm_eq(1,1) = norm_eq(1,1) + 1 bvec(1) = bvec(1) + gain_L2o(i,j) end if end if end do end do **** Now get the code-tracking L2 scaling if( norm_eq(1,1).ne.0.d0 ) then L2_code_scale = bvec(1)/norm_eq(1,1) else L2_code_scale = 1.d0 end if **** Output value write(*,230) L2_code_scale 230 format('* L2 Code Scale: ',f12.4) **** Now scale the remaining L2 code tracking gains. do i = 1, num_epochs do j = 1, num_sat if( kbit(flags(i,j),1) .and. * ! SNR good and code tracking . kbit(flags(i,j),3) .and. . num_L2_code(1).lt.num_L2_code(2) ) then gain_L2o(i,j) = gain_L2o(i,j)/L2_code_scale end if if( kbit(flags(i,j),1) .and. * ! SNR good and X-correlation . .not.kbit(flags(i,j),3) .and. . num_L2_code(1).ge.num_L2_code(2) ) then gain_L2o(i,j) = gain_L2o(i,j)/L2_code_scale end if end do end do **** Now write out the Gain curves as a function of elevation angle if( kbit(out_opts,1) ) then write(*,250) 250 format('+ Gain Curves for L1 and L2 GAIN',/, . '+ Elev SNR L1 Gain L2 SNR L2 GAIN') do elv = 90,1,-1 se = sin(elv/rad_to_deg) PG_L1 = 0 PG_L2 = 0 do k = 1, num_poly_L1 PG_L1 = PG_L1 + gain_poly_L1(k)*se**(k-1) end do do k = 1, num_poly_L2 PG_L2 = PG_L2 + gain_poly_L2(k)*se**(k-1) end do write(*,260) elv, PG_L1, PG_L2, PG_L2*PG_L1/100 260 format('+',F10.2,3F12.4,' GAIN') end do endif **** See if we should write out Normalized SNR values. if( kbit(out_opts,2) ) then write(*,300) 300 format('* Normalized SNR Observed values OBSG',/, . '* Date h m s sod PRN ', . ' Azimuth Elev L1 NSNR L2 NSNR C OBSG') do i = 1, num_epochs call mjd_to_ymdhms( epoch(i)+1.d-8, date, sectag) sod = date(4)*3600.d0 + date(5)*60.d0 + sectag do j = 1, num_sat if( kbit(flags(i,j),1) ) then code_type = 1 if( kbit(flags(i,j),3) ) code_type = 0 write(*,320) date, sectag, sod, prn_list(j), . az(i,j)*180/pi, el(i,j)*180/pi, . gain_L1o(i,j), gain_L2o(i,j), code_type 320 format(1x,i4,4i3,f6.2,1x,f8.2,' PRN ',i2.2,2f9.4, . 2f10.5,1x,i2, ' OBSG') end if end do end do end if * Thats all return end CTITLE CLEAR_NORM subroutine clear_norm( ndim ) *** Routine to clear the normal equations include 'spcsnr.h' * PASSED variables * ndim - Dimension to cleared integer*4 ndim * LOCAL variables * i,j - Loop counters integer*4 i,j **** Loop over the matrix do i = 1, ndim bvec(i) = 0.0d0 do j = 1, ndim norm_eq(i,j) = 0.0d0 end do end do **** Thats all return end CTITLE INC_NORM subroutine inc_norm( omc, apart, ndim ) * Rouitne to increment the normal equations include 'spcsnr.h' * PASSED variables * ndim - number of parameters being estimated integer*4 ndim * omc - Observed values to fit * apart(ndim) - Partial derivatives real*8 omc, apart(ndim) * LOCAL variables * i,j - Loop counters integer*4 i,j **** Loop over the system of equations do i = 1, ndim bvec(i) = bvec(i) + apart(i)*omc do j = 1, ndim norm_eq(i,j) = norm_eq(i,j) + apart(i)*apart(j) end do end do **** Thats all return end CTITLE ESTPHS subroutine estphs * This routine uses the smoothed SNR values to obtain estimates of * phase errors. The basic model is to look for increasing or decreasing * gain and to use these to get the estimates of the mean signal level * and noise amplitude. Values in-between the extremes are estimated from * the SNR value at these times. include 'spcsnr.h' include 'const_param.h' * LOCAL VARIABLES * i,j,k - Loop counters * dir - +1 if SNR is increasing, and -1 if decreasing * ep1 - Starting epoch of current sequence integer*4 i,j,k, dir, ep1 * kbit - Check status of bits * finished - Set true when satellite is completed * OK - Set true if there is good data +-2 epochs of where we * - are looking. logical kbit, finished, OK * g1, g2 - First and last gain values in current sequence * go _ Observed amplitude at epoch * gm, amp - Mean gain and the amplitude of the noise. * dphs - Correction to the phase in cycles. * cosdphs - Cosine of the phase error * deldt - Sign of elevation angle derivative (+1 when rising) real*8 g1, g2, go, gm, amp, dphs, cosdphs, deldt ***** Loop over each satellite at L1 first do j = 1, num_sat * Work our way up the first SNR value i = 1 do while ( .not.kbit(flags(i,j),1) .and. i.lt.num_epochs ) i = i + 1 end do if( gain_L1s(i+1,j).gt. gain_L1s(i,j) ) then dir = +1 else dir = -1 end if finished = .true. if( i.lt.num_epochs ) finished = .false. ep1 = i g1 = gain_L1s(ep1,j) do while ( .not.finished .and. i.lt.num_epochs ) **** See if SNR is still going in the same direction i = i + 1 if( i.ge.num_epochs ) finished = .true. ***** Check to see if we have data near this epoch OK = .false. do k = max(1,i-2), min(num_epochs,i+2) if( kbit(flags(k,j),1) ) OK = .true. end do *** Logic not quite correct. (Misses the end of the sequence if * gain does not change direction). if( dir*(gain_L1s(i+1,j)-gain_L1s(i,j)).lt.0 . .and. OK) then * The gain has changed direction. Get the mean signal * and noise amplitude and compute phase errors. g2 = gain_L1s(i,j) gm = (g1 + g2 )/2 amp = abs(g1-g2)/2 * Now loop correction the phase. (Loop to point before * i. We will then start next scan at i. do k = ep1, i-1 go = gain_L1s(k,j) * * get the sign of the elevation angle time derivative deldt = -1.d0 if( el(k+1,j).gt.el(k,j) ) deldt = 1.d0 * Use cosine rule to get the change in phase cosdphs = (go**2+gm**2-amp**2)/(2*go*gm) if( cosdphs.gt.1.d0 ) cosdphs = 1.d0 if( cosdphs.lt.-1.d0 ) cosdphs = -1.d0 dphs = acos(cosdphs)/(2*pi) phs_L1a(k,j) = phs_L1a(k,j) - dir*dphs*deldt end do ep1 = i g1 = gain_L1s(ep1,j) dir = -dir end if **** If value is not OK, then we have hit the end of sequence * so search for next start up vallue if( .not.OK ) then do while ( .not.kbit(flags(i,j),1) . .and. i.lt.num_epochs ) i = i + 1 end do if( gain_L1s(i+1,j).gt. gain_L1s(i,j) ) then dir = +1 else dir = -1 end if ep1 = i end if end do end do ***** Loop over each satellite at L2 Now do j = 1, num_sat * Work our way up the first SNR value i = 1 do while ( .not.kbit(flags(i,j),1) .and. i.lt.num_epochs ) i = i + 1 end do if( gain_L2s(i+1,j).gt. gain_L2s(i,j) ) then dir = +1 else dir = -1 end if finished = .true. if( i.lt.num_epochs ) finished = .false. ep1 = i g1 = gain_L2s(ep1,j) do while ( .not.finished .and. i.lt.num_epochs ) **** See if SNR is still going in the same direction i = i + 1 if( i.ge.num_epochs ) finished = .true. ***** Check to see if we have data near this epoch OK = .false. do k = max(1,i-2), min(num_epochs,i+2) if( kbit(flags(k,j),1) ) OK = .true. end do *** Logic not quite correct. (Misses the end of the sequence if * gain does not change direction). if( dir*(gain_L2s(i+1,j)-gain_L2s(i,j)).lt.0 . .and. OK) then * The gain has changed direction. Get the mean signal * and noise amplitude and compute phase errors. g2 = gain_L2s(i,j) gm = (g1 + g2 )/2 amp = abs(g1-g2)/2 * Now loop correction the phase. (Loop to point before * i. We will then start next scan at i. do k = ep1, i-1 go = gain_L2s(k,j) * * get the sign of the elevation angle time derivative deldt = -1.d0 if( el(k+1,j).gt.el(k,j) ) deldt = 1.d0 * Use cosine rule to get the change in phase cosdphs = (go**2+gm**2-amp**2)/(2*go*gm) if( cosdphs.gt.1.d0 ) cosdphs = 1.d0 if( cosdphs.lt.-1.d0 ) cosdphs = -1.d0 dphs = acos(cosdphs)/(2*pi) phs_L2a(k,j) = phs_L2a(k,j) - dir*dphs*deldt end do ep1 = i g1 = gain_L2s(ep1,j) dir = -dir end if **** If value is not OK, then we have hit the end of sequence * so search for next start up vallue if( .not.OK ) then do while ( .not.kbit(flags(i,j),1) . .and. i.lt.num_epochs ) i = i + 1 end do if( gain_L2s(i+1,j).gt. gain_L2s(i,j) ) then dir = +1 else dir = -1 end if ep1 = i end if end do end do **** Thats all return end CTITLE EXT_DATA subroutine ext_data( pn, more_data, split_seq , band ) * Routine to pull a sequence of data from the SNR records include 'spcsnr.h' * PASSED variables * pn - PRN number integer*4 pn * more_data - Logical to indicate more data * split_Seq - Logical to show sequence has been split because it * is long logical more_data, split_seq * band - Either L1 or L2 to indicate which one we are extracting character*(*) band * LOCAL variables * i,j - Loop counters * missed - Number of missed entries in a sequence * edbit - Indicates which bit to check for good data (4-L1, 5-L2) integer*4 i,j, missed, edbit * dg_L1, dg_L2 - change in gain over gap so that we can linearly * interpolate. * dedt - Rate of change of elevation angle. Sequence stopped when * this changes sign. * init_dedt - Angle rate of change when sequence started. * prev_el - Elevation at previous epoch. real*8 dg_L1, dg_L2, dedt, init_dedt, prev_el * start, tail - logical indicating that start and tail have been * found logical start, tail, kbit ***** Set up if( band(1:2).eq. 'L1' ) then edbit = 4 else edbit = 5 end if more_data = .true. **** Search for the start of a sequence i = start_seq start = .false. do while ( .not.start) i = i + 1 if( i.ge.num_epochs ) then more_data = .false. start = .true. end if if( kbit(flags(i,pn),edbit) ) then * Good data found start_seq = i start = .true. init_dedt = 0.d0 prev_el = el(i,pn) end if end do * If there is more data then find the end to the sequence if( more_data ) then tail = .false. missed = 0 do while ( .not.tail ) i = i + 1 * if we are about the pass the end of the data * then stop. if ( i.ge.num_epochs ) then tail = .true. end if * we allow some missing data. Keep track of how many if( .not.kbit(flags(i,pn),edbit) ) then missed = missed + 1 else end_seq = i dedt = el(i,pn)-prev_el if( init_dedt.ne.0.d0 ) then * Sign of rate of change of elevation has * changed, so end the sequence if( dedt*init_dedt.lt. 0 ) then tail = .true. end if else init_dedt = dedt end if * Save the elevation angle prev_el = el(i,pn) * Interpolate any missed data if( missed.gt.0 ) then dg_L1= gain_L1o(i,pn) - gain_L1o(i-missed-1,pn) dg_L2= gain_L2o(i,pn) - gain_L2o(i-missed-1,pn) do j = 1, missed if( band(1:2).eq.'L1' ) then gain_L1o(i-j,pn) = gain_L1o(i-missed-1,pn) + . (dg_L1/(missed+1))*j else gain_L2o(i-j,pn) = gain_L2o(i-missed-1,pn) + . (dg_L2/(missed+1))*j end if end do missed = 0 end if end if * We have missed too many epochs, so stop this sequence if( missed.gt.10 ) then tail = .true. init_dedt = 0.d0 end if end do ***** Now copy the gain values over num_seq = end_seq - start_seq + 1 if( num_seq.le.0 ) then end_seq = start_seq num_seq = 0 end if * Check the length of the sequence. If it is very long then * divide in half. if( num_seq.gt. 5000 .and. .not.split_seq ) then end_seq = start_seq + num_seq/2 num_seq = end_seq - start_seq + 1 split_seq = .true. else split_seq = .false. end if do j = 1, num_seq if( band.eq.'L1' ) then gain_L1q(j) = gain_L1o(start_seq+j-1,pn) else gain_L2q(j) = gain_L2o(start_seq+j-1,pn) end if end do if( kbit(out_opts,3) ) . write(*,200) prn_list(pn), band, start_seq, num_seq 200 format('* For PRN ',i2.2,' Band ',a,' Sequence start ', . i5,' with ', i4,' in sequence') end if ***** Thats all return end CTITLE COPY_FFT subroutine copy_fft( gain, seq, num_seq, max_seq ) * Routine to copy R*8 array to complex real*4 for FFT * PASSED variables * num_seq - number in sequence * max_seq - Maximum length of sequence integer*4 num_seq, max_seq * Gain - Gain sequence real*8 gain(max_seq) * seq - Real*4 complex copy real*4 seq(2,max_seq) * LOCAL integer*4 i do i = 1, num_seq seq(1,i) = gain(i) seq(2,i) = 0.0 end do **** Thats all return end CTITLE SWAP_SEQ subroutine swap_seq( seq, num_seq, max_seq) * Routine to swap the real and imaginary parts of the FFT * PASSED variables * num_seq - number in sequence * max_seq - Maximum length of sequence integer*4 num_seq, max_seq * seq - Real*4 complex copy real*4 seq(2,max_seq) * LOCAL integer*4 i real*4 tmp * Zero the mean before we start seq(1,1) = 0.0 seq(2,1) = 0.0 do i = 2, num_seq tmp = seq(1,i) seq(1,i) = seq(2,i) seq(2,i) = tmp end do **** Thats all return end CTITLE COMP_DPHS subroutine comp_dphs( gain, dphs, pn) * Routine to compute the phase changes from the FFT results include 'spcsnr.h' include 'const_param.h' * Passed variables integer*4 pn real*8 gain(max_seq), dphs(max_epochs, max_sat) * LOCAL variables integer*4 i, ep real*8 dp **** Loop over the current sequence and compute the corrections * to the phase. do i = 1, num_seq ep = start_seq + i - 1 if( abs(seq(2,i)).lt.gain(i) ) then dp = asin(seq(2,i)/gain(i))/(2*pi) else dp = 0 end if if( ep.lt.num_epochs ) then if( el(ep+1,pn).gt.el(ep,pn)) dp = -dp else if( el(ep,pn).gt.el(ep-1,pn)) dp = -dp end if dphs(ep,pn) = -dp end do ***** Thats all return end CTITLE FIND_MAX_SPEC subroutine find_max_spec( seq, num_seq, mf, flat, out_opts ) * Routine to find the highest peak in the first half * of the spectrum * PASSED * mf - Index of max in spectrum minus 1 (so that this is * mf/num_seq is the frequency. cycles/ep. * num_seq - Length of sequence * out_opts - Output options integer*4 mf, num_seq, out_opts * seq(2,num_seq) - FFT of SNR data. real*4 seq(2,num_seq) * Flat - Indicates that the spectrum is flat (largest peak * <1% of mean. logical*4 flat * LOCAL integer*4 i real*4 max_power, power logical kbit **** Initialize flat = .false. max_power = 0.d0 * Start at index 3 to miss mean, and the lowest frequency * term because sign of this term is not clear! * Now that we split the segments go back to low way lengths do i = 3, num_seq/4 + 1 power = seq(1,i)**2 + seq(2,i)**2 if( power.gt.max_power ) then mf = i - 1 max_power = power end if end do ***** See if spectrum is flat. Changed 1000 to 500 if( sqrt(max_power)*500 .lt. seq(1,1) ) then flat = .true. end if if( kbit(out_opts,3) ) .write(*,200) mf, sqrt(max_power),seq(1,1),float(num_seq)/mf, . flat 200 format('* MAX: Freq ',i4,' SNRs ',2f12.4,' Period ',f10.2, . ' Epochs, Flat? ',L1) **** Thats all return end CTITLE EST_SNR_SPEC Subroutine est_snr_spec( gain ) * Routine to estimate the coefficients of the SNR variations. * Separate coefficients are estimated for each cycle of SNR * frequency. Once computed, these coefficients are removed * from the gain. (Process is iterated). include 'spcsnr.h' include 'const_param.h' * PASSED * gain - Gain that we fit mean, cos, and sin. real*8 gain(num_seq) * LOCAL Variables integer*4 i,j, ipivot(3) real*8 apart(3), scale(3), period ***** Loop over each cycle in sequence. Compute period first period = (num_seq*1.d0)/mf do i = 1, mf call clear_norm( 3 ) do j = (i-1)*period+1,i*period apart(1) = 1.d0 apart(2) = cos(2*pi*j/period) apart(3) = sin(2*pi*j/period) call inc_norm( gain(j), apart, 3 ) end do ***** Solve the system of equations call invert_vis( norm_eq, bvec, scale, ipivot, 3, . max_poly, 1 ) ***** Now save the results num_cmp = num_cmp + 1 per_cmp(num_cmp) = period bnd_cmp(num_cmp) = i cos_cmp(num_cmp) = bvec(2) sin_cmp(num_cmp) = bvec(3) * Now remove contribution from Gain. do j = (i-1)*period+1,i*period gain(j) = gain(j) - . cos(2*pi*j/period)*bvec(2) - . sin(2*pi*j/period)*bvec(3) end do end do ***** Thats all. return end CTITLE SUM_DPHS subroutine sum_dphs( gain, phs, pn) * Routine to sum up all the components to get the phase * error. include 'spcsnr.h' include 'const_param.h' * PASSED * pn - PRN number being processed integer*4 pn * gain(max_epochs,max_sat) - Gain measurements. (These will be flat * by this time. * phs(max_epochs,max_sat) - Corrections to the phase measurements. real*8 gain(max_epochs,max_sat), phs(max_epochs,max_sat) * LOCAL * i,j - Loop counters * ep - Epoch number integer*4 i,j, ep * mean_sig - Mean signal (should be near 1) * sum_real, sum_imag - Sum of the real and imaginary components * of the signal * dp - Phase correction (cycles). Sign set by rate of change * of elevation angle * arg - Argument for sin and cos real*8 mean_sig, sum_real, sum_imag, dp, arg ***** Start, set the mean signal amplitude mean_sig = seq(1,1)/sqrt(1.d0*num_seq) * Now step over all the epochs in this sequence do j = 1, num_seq * Clear the real and imaginary parts sum_real = mean_sig sum_imag = 0.d0 ep = start_seq + j - 1 * Now see which components apply do i = 1, num_cmp if( j.ge.int((bnd_cmp(i)-1)*per_cmp(i)+1) .and. . j.le.int((bnd_cmp(i)*per_cmp(i))) ) then arg = 2*pi*j/per_cmp(i) sum_real = sum_real + cos_cmp(i)*cos(arg) + . sin_cmp(i)*sin(arg) sum_imag = sum_imag + sin_cmp(i)*cos(arg) - . cos_cmp(i)*sin(arg) end if end do ***** Now get the phase corrections dp = atan2(sum_imag,sum_real)/(2*pi) * Use the elevation angle rate to see satellite * going up or down. if( ep.lt.num_epochs ) then if( el(ep+1,pn).gt.el(ep,pn)) dp = -dp else if( el(ep,pn).gt.el(ep-1,pn)) dp = -dp end if phs(ep,pn) = dp gain(ep,pn) = sqrt(sum_real**2+sum_imag**2) end do **** Thats all return end CTITLE spcbd block data spcbd * Block data that sets the output options commands include 'spcsnr.h' data out_cmds / 'GAIN','OBSG','FFTC','DPHS','SMTH','RMSE', . 'SP07','SP08','SP09','SP10','SP11','SP12', . 'SP13','SP14','SP15','SP16','SP17','SP18', . 'SP19','SP20','SP21','SP22','SP23','SP24', . 'SP25','SP26','SP27','SP28','SP29','SP30', . 'SP31','SP32' / end