Browse Source

cleaning folder structure

pull/2/head
DfAC 8 years ago
parent
commit
9d41fb8fa3
  1. 95
      GNSSLogger/app/app.iml
  2. 2
      GNSSLogger/build.gradle
  3. 5
      GNSSLogger/local.properties
  4. 0
      ProcessGNSSMeas/CheckDataFilter.m
  5. 0
      ProcessGNSSMeas/CheckGpsEphInputs.m
  6. 0
      ProcessGNSSMeas/ClosestGpsEph.m
  7. 0
      ProcessGNSSMeas/CompareVersions.m
  8. 0
      ProcessGNSSMeas/Contents.m
  9. 0
      ProcessGNSSMeas/DayOfYear.m
  10. 0
      ProcessGNSSMeas/FlightTimeCorrection.m
  11. 0
      ProcessGNSSMeas/GetNasaHourlyEphemeris.m
  12. 0
      ProcessGNSSMeas/GnssThresholds.m
  13. 310
      ProcessGNSSMeas/Gps2Utc.m
  14. 0
      ProcessGNSSMeas/GpsAdrResiduals.m
  15. 0
      ProcessGNSSMeas/GpsConstants.m
  16. 178
      ProcessGNSSMeas/GpsEph2Dtsv.m
  17. 132
      ProcessGNSSMeas/GpsEph2Pvt.m
  18. 304
      ProcessGNSSMeas/GpsEph2Xyz.m
  19. 0
      ProcessGNSSMeas/GpsWlsPvt.m
  20. 108
      ProcessGNSSMeas/JulianDay.m
  21. 118
      ProcessGNSSMeas/Kepler.m
  22. 0
      ProcessGNSSMeas/LICENSE.txt
  23. 188
      ProcessGNSSMeas/LeapSeconds.m
  24. 124
      ProcessGNSSMeas/Lla2Ned.m
  25. 108
      ProcessGNSSMeas/Lla2Xyz.m
  26. 0
      ProcessGNSSMeas/PlotAdr.m
  27. 0
      ProcessGNSSMeas/PlotAdrResids.m
  28. 0
      ProcessGNSSMeas/PlotCno.m
  29. 0
      ProcessGNSSMeas/PlotPseudorangeRates.m
  30. 0
      ProcessGNSSMeas/PlotPseudoranges.m
  31. 0
      ProcessGNSSMeas/PlotPvt.m
  32. 0
      ProcessGNSSMeas/PlotPvtStates.m
  33. 0
      ProcessGNSSMeas/ProcessAdr.m
  34. 0
      ProcessGNSSMeas/ProcessGnssMeas.m
  35. 0
      ProcessGNSSMeas/ProcessGnssMeasScript.m
  36. 0
      ProcessGNSSMeas/README.md
  37. 0
      ProcessGNSSMeas/ReadGnssLogger.m
  38. 514
      ProcessGNSSMeas/ReadRinexNav.m
  39. 120
      ProcessGNSSMeas/RotEcef2Ned.m
  40. 0
      ProcessGNSSMeas/SetDataFilter.m
  41. 222
      ProcessGNSSMeas/Utc2Gps.m
  42. 0
      ProcessGNSSMeas/WlsPvt.m
  43. 164
      ProcessGNSSMeas/Xyz2Lla.m
  44. 3352
      ProcessGNSSMeas/demoFiles/hour1820.16n
  45. 0
      ProcessGNSSMeas/demoFiles/hour2570.16n
  46. 3478
      ProcessGNSSMeas/demoFiles/raw.csv
  47. 0
      ProcessGNSSMeas/demoFiles/workshop_trials01.txt

95
GNSSLogger/app/app.iml

@ -8,12 +8,12 @@
</facet>
<facet type="android" name="Android">
<configuration>
<option name="SELECTED_BUILD_VARIANT" value="debug" />
<option name="SELECTED_BUILD_VARIANT" value="release" />
<option name="SELECTED_TEST_ARTIFACT" value="_android_test_" />
<option name="ASSEMBLE_TASK_NAME" value="assembleDebug" />
<option name="COMPILE_JAVA_TASK_NAME" value="compileDebugSources" />
<option name="ASSEMBLE_TASK_NAME" value="assembleRelease" />
<option name="COMPILE_JAVA_TASK_NAME" value="compileReleaseSources" />
<afterSyncTasks>
<task>generateDebugSources</task>
<task>generateReleaseSources</task>
</afterSyncTasks>
<option name="ALLOW_USER_CONFIGURATION" value="false" />
<option name="MANIFEST_FILE_RELATIVE_PATH" value="/src/main/AndroidManifest.xml" />
@ -24,40 +24,33 @@
</facet>
</component>
<component name="NewModuleRootManager" LANGUAGE_LEVEL="JDK_1_7" inherit-compiler-output="false">
<output url="file://$MODULE_DIR$/build/intermediates/classes/debug" />
<output-test url="file://$MODULE_DIR$/build/intermediates/classes/test/debug" />
<output url="file://$MODULE_DIR$/build/intermediates/classes/release" />
<output-test url="file://$MODULE_DIR$/build/intermediates/classes/test/release" />
<exclude-output />
<content url="file://$MODULE_DIR$">
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/r/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/aidl/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/buildConfig/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/rs/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/apt/debug" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/rs/debug" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/resValues/debug" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/r/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/aidl/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/buildConfig/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/rs/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/apt/androidTest/debug" isTestSource="true" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/rs/androidTest/debug" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/resValues/androidTest/debug" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/res" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/resources" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/assets" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/aidl" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/java" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/jni" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/rs" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/debug/shaders" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/res" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/resources" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/assets" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/aidl" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/java" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/jni" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/rs" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testDebug/shaders" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/r/release" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/aidl/release" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/buildConfig/release" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/rs/release" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/source/apt/release" isTestSource="false" generated="true" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/rs/release" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/build/generated/res/resValues/release" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/release/res" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/release/resources" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/release/assets" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/release/aidl" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/release/java" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/release/jni" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/release/rs" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/release/shaders" isTestSource="false" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/res" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/resources" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/assets" type="java-test-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/aidl" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/java" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/jni" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/rs" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/testRelease/shaders" isTestSource="true" />
<sourceFolder url="file://$MODULE_DIR$/src/main/res" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/main/resources" type="java-resource" />
<sourceFolder url="file://$MODULE_DIR$/src/main/assets" type="java-resource" />
@ -84,13 +77,9 @@
<sourceFolder url="file://$MODULE_DIR$/src/androidTest/shaders" isTestSource="true" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/assets" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/blame" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/builds" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/classes" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/dependency-cache" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support.test.espresso/espresso-core/2.2.2/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support.test.espresso/espresso-idling-resource/2.2.2/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support.test/exposed-instrumentation-api-publish/0.5/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support.test/rules/0.5/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support.test/runner/0.5/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support/animated-vector-drawable/24.2.1/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support/appcompat-v7/24.2.1/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support/design/24.2.1/jars" />
@ -104,11 +93,18 @@
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support/support-v4/24.2.1/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/exploded-aar/com.android.support/support-vector-drawable/24.2.1/jars" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental-classes" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental-runtime-classes" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental-safeguard" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/incremental-verifier" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/instant-run-resources" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/instant-run-support" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/jniLibs" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/manifests" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/pre-dexed" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/reload-dex" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/res" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/restart-dex" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/rs" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/shaders" />
<excludeFolder url="file://$MODULE_DIR$/build/intermediates/symbols" />
@ -120,29 +116,18 @@
<orderEntry type="sourceFolder" forTests="false" />
<orderEntry type="library" exported="" name="support-v4-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-compat-24.2.1" level="project" />
<orderEntry type="library" exported="" name="animated-vector-drawable-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-fragment-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="runner-0.5" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="espresso-idling-resource-2.2.2" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-library-1.3" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-integration-1.3" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="jsr305-2.0.1" level="project" />
<orderEntry type="library" exported="" name="design-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-media-compat-24.2.1" level="project" />
<orderEntry type="library" exported="" name="animated-vector-drawable-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-v13-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="espresso-core-2.2.2" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="exposed-instrumentation-api-publish-0.5" level="project" />
<orderEntry type="library" exported="" name="support-fragment-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-core-ui-24.2.1" level="project" />
<orderEntry type="library" exported="" name="recyclerview-v7-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="rules-0.5" level="project" />
<orderEntry type="library" exported="" name="appcompat-v7-24.2.1" level="project" />
<orderEntry type="library" exported="" name="support-vector-drawable-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javax.annotation-api-1.2" level="project" />
<orderEntry type="library" exported="" name="support-core-utils-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javax.inject-1" level="project" />
<orderEntry type="library" exported="" name="support-annotations-24.2.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="javawriter-2.1.1" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="hamcrest-core-1.3" level="project" />
<orderEntry type="library" exported="" scope="TEST" name="junit-4.12" level="project" />
<orderEntry type="library" exported="" name="design-24.2.1" level="project" />
</component>
</module>

2
GNSSLogger/build.gradle

@ -5,7 +5,7 @@ buildscript {
jcenter()
}
dependencies {
classpath 'com.android.tools.build:gradle:2.2.0'
classpath 'com.android.tools.build:gradle:2.2.3'
// NOTE: Do not place your application dependencies here; they belong
// in the individual module build.gradle files

5
GNSSLogger/local.properties

@ -1,10 +1,11 @@
## This file is automatically generated by Android Studio.
# Do not modify this file -- YOUR CHANGES WILL BE ERASED!
#
# This file should *NOT* be checked into Version Control Systems,
# This file must *NOT* be checked into Version Control Systems,
# as it contains information specific to your local configuration.
#
# Location of the SDK. This is only used by Gradle.
# For customization when using a Version Control System, please read the
# header note.
sdk.dir=/media/build/master/prebuilts/fullsdk/linux
#Wed Mar 08 14:07:24 GMT 2017
sdk.dir=C\:\\Utils\\Android\\sdk

0
opensource/CheckDataFilter.m → ProcessGNSSMeas/CheckDataFilter.m

0
opensource/CheckGpsEphInputs.m → ProcessGNSSMeas/CheckGpsEphInputs.m

0
opensource/ClosestGpsEph.m → ProcessGNSSMeas/ClosestGpsEph.m

0
opensource/CompareVersions.m → ProcessGNSSMeas/CompareVersions.m

0
opensource/Contents.m → ProcessGNSSMeas/Contents.m

0
opensource/DayOfYear.m → ProcessGNSSMeas/DayOfYear.m

0
opensource/FlightTimeCorrection.m → ProcessGNSSMeas/FlightTimeCorrection.m

0
opensource/GetNasaHourlyEphemeris.m → ProcessGNSSMeas/GetNasaHourlyEphemeris.m

0
opensource/GnssThresholds.m → ProcessGNSSMeas/GnssThresholds.m

310
opensource/Gps2Utc.m → ProcessGNSSMeas/Gps2Utc.m

@ -1,155 +1,155 @@
function [utcTime] = Gps2Utc(gpsTime,fctSeconds)
% [utcTime] = Gps2Utc(gpsTime,[fctSeconds])
% Convert GPS time (week & seconds), or Full Cycle Time (seconds) to UTC
%
% Input: gpsTime, [mx2] matrix [gpsWeek, gpsSeconds],
% fctSeconds, [optional] Full Cycle Time (seconds)
%
% Outputs: utcTime, [mx6] matrix = [year,month,day,hours,minutes,seconds]
%
% If fctSeconds is provided, gpsTime is ignored
%
% Valid range of inputs:
% gps times corresponding to 1980/6/1 <= time < 2100/1/1
% i.e. [0,0] <= gpsTime < [6260, 432000]
% 0 <= fctSeconds < 3786480000
%
% See also: Utc2Gps
% Other functions needed: LeapSeconds.m
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Algorithm for handling leap seconds:
% When a leap second happens, utc time stands still for one second, so
% gps seconds get further ahead, so we subtract leapSecs to move from gps time
%
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
% 2) look up leap seconds for time: ls = LeapSeconds(time);
% This is (usually) the correct leap second value. Unless:
% If (utcTime=time) and (utcTime=time+ls) straddle a leap second
% then we need to add 1 to ls
% So, after step 2) ...
% 3) convert gpsTime-ls to timeMLs
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
%% Check inputs
if nargin<2 && size(gpsTime,2)~=2
error('gpsTime must have two columns')
end
if nargin<2
fctSeconds = gpsTime*[GpsConstants.WEEKSEC; 1];
end
%fct at 2100/1/1 00:00:00, not counting leap seconds:
fct2100 = [6260, 432000]*[GpsConstants.WEEKSEC; 1];
if any(fctSeconds<0) || any (fctSeconds >= fct2100)
error('gpsTime must be in this range: [0,0] <= gpsTime < [6260, 432000]');
end
%% Finished checks
%from now on work with fct
%% Apply algorithm for handling leaps seconds
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
time = Fct2Ymdhms(fctSeconds);
% 2) look up leap seconds for time: ls = LeapSeconds(time);
ls = LeapSeconds(time);
% 3) convert gpsTime-ls to timeMLs
timeMLs = Fct2Ymdhms(fctSeconds-ls);
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
if all(ls1==ls)
utcTime = timeMLs;
else
utcTime = Fct2Ymdhms(fctSeconds-ls1);
end
% NOTE:
% Gps2Utc.m doesn't produce 23:59:60, at a leap second.
% Instead, as the leap second occurs, the Gps2Utc.m sequence of
% UTC hh:mm:ss is 23:59:59, 00:00:00, 00:00:00
% and we keep it like that for code simplicity.
% Here are a sequence of UTC and GPS times around a leap second:
% formalUtcTimes = [1981 12 31 23 59 59; 1981 12 31 23 59 60; 1982 1 1 0 0 0];
% gpsTimes = [103 431999; 103 432000; 103 432001];
% >> Gps2Utc(gpsTimes)
% ans =
% 1981 12 31 23 59 59
% 1982 1 1 0 0 0
% 1982 1 1 0 0 0
% If you want to change this you could check LeapSeconds.m to see if you're
% exactly on the addition of a leap second, and then change the UTC format
% to include the '60' seconds
end %end of function Gps2Utc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function time = Fct2Ymdhms(fctSeconds)
%Utility function for Gps2Utc
%Convert GPS full cycle time to [yyyy,mm,dd,hh,mm,ss.s] format
HOURSEC = 3600; MINSEC = 60;
monthDays=[31,28,31,30,31,30,31,31,30,31,30,31];%days each month (not leap year)
m=length(fctSeconds);
days = floor(fctSeconds / GpsConstants.DAYSEC) + 6;%days since 1980/1/1
years=zeros(m,1)+1980;
%decrement days by a year at a time, until we have calculated the year:
leap=ones(m,1); %1980 was a leap year
while (any(days > (leap+365)))
I = find(days > (leap+365) );
days(I) = days(I) - (leap(I) + 365);
years(I)=years(I)+1;
leap(I) = (rem(years(I),4) == 0); % leap = 1 on a leap year, 0 otherwise
% This works from 1901 till 2099, 2100 isn't a leap year (2000 is).
% Calculate the year, ie time(1)
end,
time=zeros(m,6);%initialize matrix
time(:,1)=years;
%decrement days by a month at a time, until we have calculated the month
% Calculate the month, ie time(:,2)
% Loop through m:
for i=1:m
month = 1;
if (rem(years(i),4) == 0) %This works from 1901 till 2099
monthDays(2)=29; %Make February have 29 days in a leap year
else
monthDays(2)=28;
end
while (days(i) > monthDays(month))
days(i) = days(i)-monthDays(month);
month = month+1;
end
time(i,2)=month;
end
time(:,3) = days;
sinceMidnightSeconds = rem(fctSeconds, GpsConstants.DAYSEC);
time(:,4) = fix(sinceMidnightSeconds/HOURSEC);
lastHourSeconds = rem(sinceMidnightSeconds, HOURSEC);
time(:,5) = fix(lastHourSeconds/MINSEC);
time(:,6) = rem(lastHourSeconds,MINSEC);
end %end of function Fct2Ymdhms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [utcTime] = Gps2Utc(gpsTime,fctSeconds)
% [utcTime] = Gps2Utc(gpsTime,[fctSeconds])
% Convert GPS time (week & seconds), or Full Cycle Time (seconds) to UTC
%
% Input: gpsTime, [mx2] matrix [gpsWeek, gpsSeconds],
% fctSeconds, [optional] Full Cycle Time (seconds)
%
% Outputs: utcTime, [mx6] matrix = [year,month,day,hours,minutes,seconds]
%
% If fctSeconds is provided, gpsTime is ignored
%
% Valid range of inputs:
% gps times corresponding to 1980/6/1 <= time < 2100/1/1
% i.e. [0,0] <= gpsTime < [6260, 432000]
% 0 <= fctSeconds < 3786480000
%
% See also: Utc2Gps
% Other functions needed: LeapSeconds.m
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Algorithm for handling leap seconds:
% When a leap second happens, utc time stands still for one second, so
% gps seconds get further ahead, so we subtract leapSecs to move from gps time
%
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
% 2) look up leap seconds for time: ls = LeapSeconds(time);
% This is (usually) the correct leap second value. Unless:
% If (utcTime=time) and (utcTime=time+ls) straddle a leap second
% then we need to add 1 to ls
% So, after step 2) ...
% 3) convert gpsTime-ls to timeMLs
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
%% Check inputs
if nargin<2 && size(gpsTime,2)~=2
error('gpsTime must have two columns')
end
if nargin<2
fctSeconds = gpsTime*[GpsConstants.WEEKSEC; 1];
end
%fct at 2100/1/1 00:00:00, not counting leap seconds:
fct2100 = [6260, 432000]*[GpsConstants.WEEKSEC; 1];
if any(fctSeconds<0) || any (fctSeconds >= fct2100)
error('gpsTime must be in this range: [0,0] <= gpsTime < [6260, 432000]');
end
%% Finished checks
%from now on work with fct
%% Apply algorithm for handling leaps seconds
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
time = Fct2Ymdhms(fctSeconds);
% 2) look up leap seconds for time: ls = LeapSeconds(time);
ls = LeapSeconds(time);
% 3) convert gpsTime-ls to timeMLs
timeMLs = Fct2Ymdhms(fctSeconds-ls);
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
if all(ls1==ls)
utcTime = timeMLs;
else
utcTime = Fct2Ymdhms(fctSeconds-ls1);
end
% NOTE:
% Gps2Utc.m doesn't produce 23:59:60, at a leap second.
% Instead, as the leap second occurs, the Gps2Utc.m sequence of
% UTC hh:mm:ss is 23:59:59, 00:00:00, 00:00:00
% and we keep it like that for code simplicity.
% Here are a sequence of UTC and GPS times around a leap second:
% formalUtcTimes = [1981 12 31 23 59 59; 1981 12 31 23 59 60; 1982 1 1 0 0 0];
% gpsTimes = [103 431999; 103 432000; 103 432001];
% >> Gps2Utc(gpsTimes)
% ans =
% 1981 12 31 23 59 59
% 1982 1 1 0 0 0
% 1982 1 1 0 0 0
% If you want to change this you could check LeapSeconds.m to see if you're
% exactly on the addition of a leap second, and then change the UTC format
% to include the '60' seconds
end %end of function Gps2Utc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function time = Fct2Ymdhms(fctSeconds)
%Utility function for Gps2Utc
%Convert GPS full cycle time to [yyyy,mm,dd,hh,mm,ss.s] format
HOURSEC = 3600; MINSEC = 60;
monthDays=[31,28,31,30,31,30,31,31,30,31,30,31];%days each month (not leap year)
m=length(fctSeconds);
days = floor(fctSeconds / GpsConstants.DAYSEC) + 6;%days since 1980/1/1
years=zeros(m,1)+1980;
%decrement days by a year at a time, until we have calculated the year:
leap=ones(m,1); %1980 was a leap year
while (any(days > (leap+365)))
I = find(days > (leap+365) );
days(I) = days(I) - (leap(I) + 365);
years(I)=years(I)+1;
leap(I) = (rem(years(I),4) == 0); % leap = 1 on a leap year, 0 otherwise
% This works from 1901 till 2099, 2100 isn't a leap year (2000 is).
% Calculate the year, ie time(1)
end,
time=zeros(m,6);%initialize matrix
time(:,1)=years;
%decrement days by a month at a time, until we have calculated the month
% Calculate the month, ie time(:,2)
% Loop through m:
for i=1:m
month = 1;
if (rem(years(i),4) == 0) %This works from 1901 till 2099
monthDays(2)=29; %Make February have 29 days in a leap year
else
monthDays(2)=28;
end
while (days(i) > monthDays(month))
days(i) = days(i)-monthDays(month);
month = month+1;
end
time(i,2)=month;
end
time(:,3) = days;
sinceMidnightSeconds = rem(fctSeconds, GpsConstants.DAYSEC);
time(:,4) = fix(sinceMidnightSeconds/HOURSEC);
lastHourSeconds = rem(sinceMidnightSeconds, HOURSEC);
time(:,5) = fix(lastHourSeconds/MINSEC);
time(:,6) = rem(lastHourSeconds,MINSEC);
end %end of function Fct2Ymdhms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/GpsAdrResiduals.m → ProcessGNSSMeas/GpsAdrResiduals.m

0
opensource/GpsConstants.m → ProcessGNSSMeas/GpsConstants.m

178
opensource/GpsEph2Dtsv.m → ProcessGNSSMeas/GpsEph2Dtsv.m

@ -1,89 +1,89 @@
function [dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%[dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%
% Calculate satellite clock bias
% Inputs:
% gpsEph = vector of GPS ephemeris structures defined in ReadRinexNav
% tS GPS time of week (secs) at time of transmission, calculate this as:
% tS = trx - PR/c,
% where trx = time of reception (receiver clock)
% PR = pseudorange
%
% tS may be a vector, gpsEph is a structure or vector of structures
% If gpsEph is a vector then t must have the same number of rows as gpsEph
% tS(i) is interpreted as being associated with gpsEph(i,:)
% If gpsEph is a single structure then tS may be any length
%
% Outputs:
% dtsvS = sat clock bias (seconds). GPS time = satellite time - dtsvS
% length(dtsvS) = length(tS)
%
% See IS-GPS-200 and RTCM Paper 99-92 SC104-88 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%% Check size of inputs
if min(size(tS))>1
error('tS must be a vector or a scalar, not a matrix')
end
tS=tS(:)';%make t a row vector, to match [gpsEph.*] vectors below
pt=length(tS);
[p]=length(gpsEph);
if (p>1 && pt~=p),
error('If gpsEph is a vector tS must be a vector with #rows = length(gpsEph),\n')
end
%%
%%Extract the necessary variables from gpsEph
TGD = [gpsEph.TGD];
Toc = [gpsEph.Toc];
af2 = [gpsEph.af2];
af1 = [gpsEph.af1];
af0 = [gpsEph.af0];
Delta_n = [gpsEph.Delta_n];
M0 = [gpsEph.M0];
e = [gpsEph.e];
Asqrt = [gpsEph.Asqrt];
Toe = [gpsEph.Toe];
%%
%%Calculate dependent variables ------------------------------------------------
tk = tS - Toe; %time since time of applicability
I = find(tk > 302400.0);
if any(I), tk(I) = tk(I)-GpsConstants.WEEKSEC; end,
I = find(tk < -302400.0);
if (I), tk(I) = tk(I)+GpsConstants.WEEKSEC; end,
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%%
%%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1) ----------------
dt = tS - Toc;
I = find(dt > 302400.0);
if any(I), dt(I) = dt(I)-GpsConstants.WEEKSEC; end,
I = find(dt < -302400.0);
if (I), dt(I) = dt(I)+GpsConstants.WEEKSEC; end,
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin(Ek) -TGD;
end % end of function GpsEph2Dtsv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%[dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%
% Calculate satellite clock bias
% Inputs:
% gpsEph = vector of GPS ephemeris structures defined in ReadRinexNav
% tS GPS time of week (secs) at time of transmission, calculate this as:
% tS = trx - PR/c,
% where trx = time of reception (receiver clock)
% PR = pseudorange
%
% tS may be a vector, gpsEph is a structure or vector of structures
% If gpsEph is a vector then t must have the same number of rows as gpsEph
% tS(i) is interpreted as being associated with gpsEph(i,:)
% If gpsEph is a single structure then tS may be any length
%
% Outputs:
% dtsvS = sat clock bias (seconds). GPS time = satellite time - dtsvS
% length(dtsvS) = length(tS)
%
% See IS-GPS-200 and RTCM Paper 99-92 SC104-88 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%% Check size of inputs
if min(size(tS))>1
error('tS must be a vector or a scalar, not a matrix')
end
tS=tS(:)';%make t a row vector, to match [gpsEph.*] vectors below
pt=length(tS);
[p]=length(gpsEph);
if (p>1 && pt~=p),
error('If gpsEph is a vector tS must be a vector with #rows = length(gpsEph),\n')
end
%%
%%Extract the necessary variables from gpsEph
TGD = [gpsEph.TGD];
Toc = [gpsEph.Toc];
af2 = [gpsEph.af2];
af1 = [gpsEph.af1];
af0 = [gpsEph.af0];
Delta_n = [gpsEph.Delta_n];
M0 = [gpsEph.M0];
e = [gpsEph.e];
Asqrt = [gpsEph.Asqrt];
Toe = [gpsEph.Toe];
%%
%%Calculate dependent variables ------------------------------------------------
tk = tS - Toe; %time since time of applicability
I = find(tk > 302400.0);
if any(I), tk(I) = tk(I)-GpsConstants.WEEKSEC; end,
I = find(tk < -302400.0);
if (I), tk(I) = tk(I)+GpsConstants.WEEKSEC; end,
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%%
%%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1) ----------------
dt = tS - Toc;
I = find(dt > 302400.0);
if any(I), dt(I) = dt(I)-GpsConstants.WEEKSEC; end,
I = find(dt < -302400.0);
if (I), dt(I) = dt(I)+GpsConstants.WEEKSEC; end,
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin(Ek) -TGD;
end % end of function GpsEph2Dtsv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

132
opensource/GpsEph2Pvt.m → ProcessGNSSMeas/GpsEph2Pvt.m

@ -1,66 +1,66 @@
function [xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%[xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%
% Calculate sv coordinates, in ECEF frame, sv clock bias, and sv velocity
% Inputs:
% gpsEph: vector of ephemeris structures, as defined in ReadRinexEph.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using eph2dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% vMps = [i,j,k] matrix of satellite velocity (ecef m/s)
% dtsvSDot = vector of satellite clock error rate (seconds/second)
% The row dimension of xM, dtsvS, vMps, dtsvSDot = length(gpsEph)
%
% xM and vMps are the satellite positions and velocities
% at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xM & vMps in ecef coord at time of reception
%
% functions called: GpsEph2Xyz.m
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
vMps = []; dtsvSDot = [];
[xM,dtsvS]=GpsEph2Xyz(gpsEph,gpsTime);
if isempty(xM)
return
end
%compute velocity from delta position & dtsvS at (t+0.5) - (t-0.5)
%This is better than differentiating, because both the orbital and relativity
%terms have nonlinearities that are not easily differentiable
t1 = [gpsTime(:,1), gpsTime(:,2)+0.5]; %time + 0.5 seconds
[xPlus,dtsvPlus] = GpsEph2Xyz(gpsEph,t1);
t1 = [gpsTime(:,1), gpsTime(:,2)-0.5]; %time - 0.5 seconds
[xMinus,dtsvMinus]= GpsEph2Xyz(gpsEph,t1);
vMps = xPlus - xMinus;
dtsvSDot = dtsvPlus - dtsvMinus;
end %end of function GpsEph2Pvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%[xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%
% Calculate sv coordinates, in ECEF frame, sv clock bias, and sv velocity
% Inputs:
% gpsEph: vector of ephemeris structures, as defined in ReadRinexEph.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using eph2dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% vMps = [i,j,k] matrix of satellite velocity (ecef m/s)
% dtsvSDot = vector of satellite clock error rate (seconds/second)
% The row dimension of xM, dtsvS, vMps, dtsvSDot = length(gpsEph)
%
% xM and vMps are the satellite positions and velocities
% at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xM & vMps in ecef coord at time of reception
%
% functions called: GpsEph2Xyz.m
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
vMps = []; dtsvSDot = [];
[xM,dtsvS]=GpsEph2Xyz(gpsEph,gpsTime);
if isempty(xM)
return
end
%compute velocity from delta position & dtsvS at (t+0.5) - (t-0.5)
%This is better than differentiating, because both the orbital and relativity
%terms have nonlinearities that are not easily differentiable
t1 = [gpsTime(:,1), gpsTime(:,2)+0.5]; %time + 0.5 seconds
[xPlus,dtsvPlus] = GpsEph2Xyz(gpsEph,t1);
t1 = [gpsTime(:,1), gpsTime(:,2)-0.5]; %time - 0.5 seconds
[xMinus,dtsvMinus]= GpsEph2Xyz(gpsEph,t1);
vMps = xPlus - xMinus;
dtsvSDot = dtsvPlus - dtsvMinus;
end %end of function GpsEph2Pvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

304
opensource/GpsEph2Xyz.m → ProcessGNSSMeas/GpsEph2Xyz.m

@ -1,152 +1,152 @@
function [xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
%[xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
% Calculate sv coordinates, in ECEF frame, at ttx = gpsTime
%
% Inputs:
% gpsEph: vector of GPS ephemeris structures, as defined in ReadRinexNav.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using Eph2Dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xyzM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% The row dimension of xyzM, dtsvS = length(gpsEph)
%
% xyzM = sat positions at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xyzM & v in ecef coords at time of reception
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
xyzM=[]; dtsvS=[];
[bOk,gpsEph,gpsWeek,ttxSec] = CheckGpsEphInputs(gpsEph,gpsTime);
if ~bOk
return
end
p=length(gpsEph);
%Now we are done checking and manipulating the inputs
% the time vectors: gpsWeek, ttxSec are the same length as gpsEph
%-------------------------------------------------------------------------------
%set fitIntervalSeconds
fitIntervalHours = [gpsEph.Fit_interval]';
%Rinex says "Zero if not known", so adjust for zeros
fitIntervalHours(fitIntervalHours == 0) = 2;
fitIntervalSeconds = fitIntervalHours*3600;
%Extract variables from gpsEph, into column vectors
%orbit variable names follow RINEX 2.1 Table A4
%clock variables af0, af1, af2 follow IS GPS 200
TGD = [gpsEph.TGD]';
Toc = [gpsEph.Toc]';
af2 = [gpsEph.af2]';
af1 = [gpsEph.af1]';
af0 = [gpsEph.af0]';
Crs = [gpsEph.Crs]';
Delta_n = [gpsEph.Delta_n]';
M0 = [gpsEph.M0]';
Cuc = [gpsEph.Cuc]';
e = [gpsEph.e]';
Cus = [gpsEph.Cus]';
Asqrt = [gpsEph.Asqrt]';
Toe = [gpsEph.Toe]';
Cic = [gpsEph.Cic]';
OMEGA = [gpsEph.OMEGA]';
Cis = [gpsEph.Cis]';
i0 = [gpsEph.i0]';
Crc = [gpsEph.Crc]';
omega = [gpsEph.omega]';
OMEGA_DOT=[gpsEph.OMEGA_DOT]';
IDOT = [gpsEph.IDOT]';
ephGpsWeek = [gpsEph.GPS_Week]';
%Calculate dependent variables -------------------------------------------------
%Time since time of applicability accounting for weeks and therefore rollovers
%subtract weeks first, to avoid precision errors:
tk =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toe);
I = find(abs(tk)>fitIntervalSeconds);
if ~isempty(I)
numTimes = length(I);
fprintf(sprintf('WARNING in GpsEph2Xyz.m, %d times outside fit interval.',...
numTimes));
end
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
h = sqrt(A.*(1-e.^2).*GpsConstants.mu);
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1)
%subtract weeks first, to avoid precision errors:
dt =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toc);
%Calculate satellite clock bias
sin_Ek=sin(Ek);
cos_Ek=cos(Ek);
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin_Ek -TGD;
%true anomaly:
vk=atan2(sqrt(1-e.^2).*sin_Ek./(1-e.*cos_Ek),(cos_Ek-e)./(1-e.*cos_Ek));
Phik=vk + omega; %Argument of latitude
sin_2Phik=sin(2*Phik);
cos_2Phik=cos(2*Phik);
% The next three terms are the second harmonic perturbations
duk = Cus.*sin_2Phik + Cuc.*cos_2Phik; %Argument of latitude correction
drk = Crc.*cos_2Phik + Crs.*sin_2Phik; %Radius Correction
dik = Cic.*cos_2Phik + Cis.*sin_2Phik; %Correction to Inclination
uk = Phik + duk; %Corrected argument of latitude
rk = A.*((1-e.^2)./(1+e.*cos(vk))) + drk; %Corrected radius
ik = i0 + IDOT.*tk + dik; %Corrected inclination
sin_uk=sin(uk);
cos_uk=cos(uk);
xkp = rk.*cos_uk; % Position in orbital plane
ykp = rk.*sin_uk; % Position in orbital plane
% Wk = corrected longitude of ascending node,
Wk = OMEGA + (OMEGA_DOT - GpsConstants.WE).*tk - GpsConstants.WE*[gpsEph.Toe]';
%for dtflight, see FlightTimeCorrection.m
sin_Wk = sin(Wk);
cos_Wk = cos(Wk);
xyzM=zeros(p,3);
% The matrix xyzM contains the ECEF coordinates of position
sin_ik=sin(ik);
cos_ik=cos(ik);
xyzM(:,1)=xkp.*cos_Wk-ykp.*cos_ik.*sin_Wk;
xyzM(:,2)=xkp.*sin_Wk+ykp.*cos_ik.*cos_Wk;
xyzM(:,3)=ykp.*sin_ik;
end %end of function GpsEph2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
%[xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
% Calculate sv coordinates, in ECEF frame, at ttx = gpsTime
%
% Inputs:
% gpsEph: vector of GPS ephemeris structures, as defined in ReadRinexNav.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using Eph2Dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xyzM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% The row dimension of xyzM, dtsvS = length(gpsEph)
%
% xyzM = sat positions at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xyzM & v in ecef coords at time of reception
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
xyzM=[]; dtsvS=[];
[bOk,gpsEph,gpsWeek,ttxSec] = CheckGpsEphInputs(gpsEph,gpsTime);
if ~bOk
return
end
p=length(gpsEph);
%Now we are done checking and manipulating the inputs
% the time vectors: gpsWeek, ttxSec are the same length as gpsEph
%-------------------------------------------------------------------------------
%set fitIntervalSeconds
fitIntervalHours = [gpsEph.Fit_interval]';
%Rinex says "Zero if not known", so adjust for zeros
fitIntervalHours(fitIntervalHours == 0) = 2;
fitIntervalSeconds = fitIntervalHours*3600;
%Extract variables from gpsEph, into column vectors
%orbit variable names follow RINEX 2.1 Table A4
%clock variables af0, af1, af2 follow IS GPS 200
TGD = [gpsEph.TGD]';
Toc = [gpsEph.Toc]';
af2 = [gpsEph.af2]';
af1 = [gpsEph.af1]';
af0 = [gpsEph.af0]';
Crs = [gpsEph.Crs]';
Delta_n = [gpsEph.Delta_n]';
M0 = [gpsEph.M0]';
Cuc = [gpsEph.Cuc]';
e = [gpsEph.e]';
Cus = [gpsEph.Cus]';
Asqrt = [gpsEph.Asqrt]';
Toe = [gpsEph.Toe]';
Cic = [gpsEph.Cic]';
OMEGA = [gpsEph.OMEGA]';
Cis = [gpsEph.Cis]';
i0 = [gpsEph.i0]';
Crc = [gpsEph.Crc]';
omega = [gpsEph.omega]';
OMEGA_DOT=[gpsEph.OMEGA_DOT]';
IDOT = [gpsEph.IDOT]';
ephGpsWeek = [gpsEph.GPS_Week]';
%Calculate dependent variables -------------------------------------------------
%Time since time of applicability accounting for weeks and therefore rollovers
%subtract weeks first, to avoid precision errors:
tk =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toe);
I = find(abs(tk)>fitIntervalSeconds);
if ~isempty(I)
numTimes = length(I);
fprintf(sprintf('WARNING in GpsEph2Xyz.m, %d times outside fit interval.',...
numTimes));
end
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
h = sqrt(A.*(1-e.^2).*GpsConstants.mu);
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1)
%subtract weeks first, to avoid precision errors:
dt =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toc);
%Calculate satellite clock bias
sin_Ek=sin(Ek);
cos_Ek=cos(Ek);
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin_Ek -TGD;
%true anomaly:
vk=atan2(sqrt(1-e.^2).*sin_Ek./(1-e.*cos_Ek),(cos_Ek-e)./(1-e.*cos_Ek));
Phik=vk + omega; %Argument of latitude
sin_2Phik=sin(2*Phik);
cos_2Phik=cos(2*Phik);
% The next three terms are the second harmonic perturbations
duk = Cus.*sin_2Phik + Cuc.*cos_2Phik; %Argument of latitude correction
drk = Crc.*cos_2Phik + Crs.*sin_2Phik; %Radius Correction
dik = Cic.*cos_2Phik + Cis.*sin_2Phik; %Correction to Inclination
uk = Phik + duk; %Corrected argument of latitude
rk = A.*((1-e.^2)./(1+e.*cos(vk))) + drk; %Corrected radius
ik = i0 + IDOT.*tk + dik; %Corrected inclination
sin_uk=sin(uk);
cos_uk=cos(uk);
xkp = rk.*cos_uk; % Position in orbital plane
ykp = rk.*sin_uk; % Position in orbital plane
% Wk = corrected longitude of ascending node,
Wk = OMEGA + (OMEGA_DOT - GpsConstants.WE).*tk - GpsConstants.WE*[gpsEph.Toe]';
%for dtflight, see FlightTimeCorrection.m
sin_Wk = sin(Wk);
cos_Wk = cos(Wk);
xyzM=zeros(p,3);
% The matrix xyzM contains the ECEF coordinates of position
sin_ik=sin(ik);
cos_ik=cos(ik);
xyzM(:,1)=xkp.*cos_Wk-ykp.*cos_ik.*sin_Wk;
xyzM(:,2)=xkp.*sin_Wk+ykp.*cos_ik.*cos_Wk;
xyzM(:,3)=ykp.*sin_ik;
end %end of function GpsEph2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/GpsWlsPvt.m → ProcessGNSSMeas/GpsWlsPvt.m

108
opensource/JulianDay.m → ProcessGNSSMeas/JulianDay.m

@ -1,54 +1,54 @@
function jDays = JulianDay(utcTime)
% jDays = JulianDay(utcTime);
%
% input: utcTime [mx6] matrix [year,month,day,hours,minutes,seconds]
%
% output: totalDays in Julian Days [mx1] vector (real number of days)
%
% Valid input range: 1900 < year < 2100
%Algorithm from Meeus, (1991) Astronomical Algorithms,
%see http://www.geoastro.de/elevaz/basics/meeus.htm for online summary
% valid range 1900/3/1 to 2100/2/28
% but we limit inputs to 1901 through 2099, because it's simpler
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
y = utcTime(:,1);
m = utcTime(:,2);
d = utcTime(:,3);
h = utcTime(:,4) + utcTime(:,5)/60 + utcTime(:,5)/3600;
%check that date is in valid range
if ( any(y<1901) || any (y>2099) )
error('utcTime(:,1) not in allowed range: 1900 < year < 2100')
end
i2 = m<=2; %index into months <=2
m(i2) = m(i2)+12;
y(i2) = y(i2)-1;
jDays = floor(365.25*y) + floor(30.6001*(m+1)) - 15 + 1720996.5 + d + h/24;
end %end of function JulianDay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function jDays = JulianDay(utcTime)
% jDays = JulianDay(utcTime);
%
% input: utcTime [mx6] matrix [year,month,day,hours,minutes,seconds]
%
% output: totalDays in Julian Days [mx1] vector (real number of days)
%
% Valid input range: 1900 < year < 2100
%Algorithm from Meeus, (1991) Astronomical Algorithms,
%see http://www.geoastro.de/elevaz/basics/meeus.htm for online summary
% valid range 1900/3/1 to 2100/2/28
% but we limit inputs to 1901 through 2099, because it's simpler
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
y = utcTime(:,1);
m = utcTime(:,2);
d = utcTime(:,3);
h = utcTime(:,4) + utcTime(:,5)/60 + utcTime(:,5)/3600;
%check that date is in valid range
if ( any(y<1901) || any (y>2099) )
error('utcTime(:,1) not in allowed range: 1900 < year < 2100')
end
i2 = m<=2; %index into months <=2
m(i2) = m(i2)+12;
y(i2) = y(i2)-1;
jDays = floor(365.25*y) + floor(30.6001*(m+1)) - 15 + 1720996.5 + d + h/24;
end %end of function JulianDay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

118
opensource/Kepler.m → ProcessGNSSMeas/Kepler.m

@ -1,59 +1,59 @@
function ek = Kepler(mk,e)
% ek = Kepler(mk,e)
% Kepler - Solves Kepler's equation for ek through iteration.
%
% Inputs: mk: mean anomaly (rad)
% e: eccentricity
% Outputs: ek: eccentric anomaly
%
% NOTE: mk and e may be a vectors of the same dimensions, or one may be a scalar
% the output is a vector of the same dimensions as the input vector
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Check inputs size
if min(size(mk))>1
error('mk must be a vector or a scalar, not a matrix')
end
if min(size(e))>1
error('e must be a vector or a scalar, not a matrix')
end
if length(mk)>1 && length(e)>1 && any(size(mk)~=size(e))
%both are vectors, they must have the same dimensions
error('If mk and e are both vectors they must have the same dimensions')
end
err = 1;
ek=mk;
iterCount = 0;
maxIterCount = 20;
while any(abs(err) > 1e-8) && iterCount < maxIterCount
err = ek - mk - (e.*sin(ek));
ek = ek - err;
iterCount = iterCount + 1;
if iterCount == maxIterCount
fprintf('Failed convergence on Kepler''s equation.\n')
end
end
end %end of function Kepler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function ek = Kepler(mk,e)
% ek = Kepler(mk,e)
% Kepler - Solves Kepler's equation for ek through iteration.
%
% Inputs: mk: mean anomaly (rad)
% e: eccentricity
% Outputs: ek: eccentric anomaly
%
% NOTE: mk and e may be a vectors of the same dimensions, or one may be a scalar
% the output is a vector of the same dimensions as the input vector
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Check inputs size
if min(size(mk))>1
error('mk must be a vector or a scalar, not a matrix')
end
if min(size(e))>1
error('e must be a vector or a scalar, not a matrix')
end
if length(mk)>1 && length(e)>1 && any(size(mk)~=size(e))
%both are vectors, they must have the same dimensions
error('If mk and e are both vectors they must have the same dimensions')
end
err = 1;
ek=mk;
iterCount = 0;
maxIterCount = 20;
while any(abs(err) > 1e-8) && iterCount < maxIterCount
err = ek - mk - (e.*sin(ek));
ek = ek - err;
iterCount = iterCount + 1;
if iterCount == maxIterCount
fprintf('Failed convergence on Kepler''s equation.\n')
end
end
end %end of function Kepler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/LICENSE.txt → ProcessGNSSMeas/LICENSE.txt

188
opensource/LeapSeconds.m → ProcessGNSSMeas/LeapSeconds.m

@ -1,94 +1,94 @@
function [leapSecs] = LeapSeconds(utcTime)
% leapSecs = LeapSeconds(utcTime)
% find the number of leap seconds since the GPS Epoch
%
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Output: leapSecs,
% leapSecs(i) = number of leap seconds between the GPS Epoch and utcTime(i,:)
%
% The function looks up the number of leap seconds from a UTC time table
%
% LATEST LEAP SECOND IN THE TABLE = 31 Dec 2016.
% On 1 Jan 2017: GPS-UTC=18s
% See IERS Bulletin C, https://hpiers.obspm.fr/iers/bul/bulc/bulletinc.dat
% and http://tycho.usno.navy.mil/leapsec.html
%
% Aren't Leap Seconds a pain? Yes, and very costly. Ban the Leap Second:
% Leap seconds occur on average once every two years.
% What would happen if we had no leap seconds:
% 1) Thousand of engineers would NOT spend several weeks each two years fixing
% or planning for leap second bugs.
% 2) About seven thousand years from now, solar time would be 1 hour behind UTC
% and we would need a 1 hour adjustment - similar to daylight savings time.
% 3) GMT (which is solar time) will lose its significance.
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
[m,n] = size(utcTime);
if n~=6
error('utcTime input must have 6 columns');
end
% UTC table contains UTC times (in the form of [year,month,day,hours,mins,secs])
% At each of these times a leap second had just occurred
utcTable = [1982 1 1 0 0 0;
1982 7 1 0 0 0;
1983 7 1 0 0 0;
1985 7 1 0 0 0;
1988 1 1 0 0 0;
1990 1 1 0 0 0;
1991 1 1 0 0 0;
1992 7 1 0 0 0;
1993 7 1 0 0 0;
1994 7 1 0 0 0;
1996 1 1 0 0 0;
1997 7 1 0 0 0;
1999 1 1 0 0 0;
2006 1 1 0 0 0;
2009 1 1 0 0 0;
2012 7 1 0 0 0;
2015 7 1 0 0 0;
2017 1 1 0 0 0
];
%when a new leap second is announced in IERS Bulletin C
%update the table with the UTC time right after the new leap second
tableJDays = JulianDay(utcTable) - GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%tableSeconds = tableJDays*GpsConstants.DAYSEC + utcTable(:,4:6)*[3600;60;1];
%NOTE: JulianDay returns a realed value number, corresponding to days and
% fractions thereof, so we multiply it by DAYSEC to get the full time in seconds
tableSeconds = tableJDays*GpsConstants.DAYSEC;
jDays = JulianDay(utcTime)- GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%timeSeconds = jDays*GpsConstants.DAYSEC + utcTime(:,4:6)*[3600;60;1];
timeSeconds = jDays*GpsConstants.DAYSEC;
% tableSeconds and timeSeconds now contain number of seconds since the GPS epoch
leapSecs=zeros(m,1);
for i=1:m
%add up the number of leap seconds that have occured by timeSeconds(i)
leapSecs(i) = sum(tableSeconds<=timeSeconds(i));
end
end %end of function LeapSeconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [leapSecs] = LeapSeconds(utcTime)
% leapSecs = LeapSeconds(utcTime)
% find the number of leap seconds since the GPS Epoch
%
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Output: leapSecs,
% leapSecs(i) = number of leap seconds between the GPS Epoch and utcTime(i,:)
%
% The function looks up the number of leap seconds from a UTC time table
%
% LATEST LEAP SECOND IN THE TABLE = 31 Dec 2016.
% On 1 Jan 2017: GPS-UTC=18s
% See IERS Bulletin C, https://hpiers.obspm.fr/iers/bul/bulc/bulletinc.dat
% and http://tycho.usno.navy.mil/leapsec.html
%
% Aren't Leap Seconds a pain? Yes, and very costly. Ban the Leap Second:
% Leap seconds occur on average once every two years.
% What would happen if we had no leap seconds:
% 1) Thousand of engineers would NOT spend several weeks each two years fixing
% or planning for leap second bugs.
% 2) About seven thousand years from now, solar time would be 1 hour behind UTC
% and we would need a 1 hour adjustment - similar to daylight savings time.
% 3) GMT (which is solar time) will lose its significance.
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
[m,n] = size(utcTime);
if n~=6
error('utcTime input must have 6 columns');
end
% UTC table contains UTC times (in the form of [year,month,day,hours,mins,secs])
% At each of these times a leap second had just occurred
utcTable = [1982 1 1 0 0 0;
1982 7 1 0 0 0;
1983 7 1 0 0 0;
1985 7 1 0 0 0;
1988 1 1 0 0 0;
1990 1 1 0 0 0;
1991 1 1 0 0 0;
1992 7 1 0 0 0;
1993 7 1 0 0 0;
1994 7 1 0 0 0;
1996 1 1 0 0 0;
1997 7 1 0 0 0;
1999 1 1 0 0 0;
2006 1 1 0 0 0;
2009 1 1 0 0 0;
2012 7 1 0 0 0;
2015 7 1 0 0 0;
2017 1 1 0 0 0
];
%when a new leap second is announced in IERS Bulletin C
%update the table with the UTC time right after the new leap second
tableJDays = JulianDay(utcTable) - GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%tableSeconds = tableJDays*GpsConstants.DAYSEC + utcTable(:,4:6)*[3600;60;1];
%NOTE: JulianDay returns a realed value number, corresponding to days and
% fractions thereof, so we multiply it by DAYSEC to get the full time in seconds
tableSeconds = tableJDays*GpsConstants.DAYSEC;
jDays = JulianDay(utcTime)- GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%timeSeconds = jDays*GpsConstants.DAYSEC + utcTime(:,4:6)*[3600;60;1];
timeSeconds = jDays*GpsConstants.DAYSEC;
% tableSeconds and timeSeconds now contain number of seconds since the GPS epoch
leapSecs=zeros(m,1);
for i=1:m
%add up the number of leap seconds that have occured by timeSeconds(i)
leapSecs(i) = sum(tableSeconds<=timeSeconds(i));
end
end %end of function LeapSeconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

124
opensource/Lla2Ned.m → ProcessGNSSMeas/Lla2Ned.m

@ -1,62 +1,62 @@
function [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% function to difference latitude, longitude and altitude
% and provide an answer in NED coordinates in meters.
%
% Inputs: lla1DegDegM: mx3 matrix, [latitude(deg),longitude(deg),altitude(m)]
% lla2DegDegM: mx3 or 1x3, [latitude(deg),longitude(deg),altitude(m)]
% Output: nedM = lla1DegDegM - lla2DegDegM in NED coords (meters),
%
% Useful rules of thumb for quick conversions:
% 1e-5 (5th decimal place) of a degree of latitude approx= 1.1 meters
% 1e-5 of a degree of longitude approx= cos(latitude) * 1.1 meters
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if nargin<2, error('Two inputs arguments needed'),end
[m1,n1]=size(lla1DegDegM);
[m2,n2]=size(lla2DegDegM);
if m2==1
lla2DegDegM=ones(m1,1)*lla2DegDegM;
else
if m1~=m2,
error('Second input must have one row or same number of rows as first input')
end
end
if n1~=3 || n2~=3, error('Both input matrices must have 3 columns'),end
[xyz1M] = Lla2Xyz(lla1DegDegM);
[xyz2M] = Lla2Xyz(lla2DegDegM);
refXyz = (xyz1M+xyz2M)/2;
[llaDegDegM] = Xyz2Lla(refXyz);
northM = zeros(m1,1);
eastM = zeros(m1,1);
for i=1:m1
Ce2n = RotEcef2Ned(llaDegDegM(i,1),llaDegDegM(i,2));
v = Ce2n*(xyz1M(i,:)-xyz2M(i,:))';
northM(i)=v(1);
eastM(i)=v(2);
end
downM = -lla1DegDegM(:,3)+lla2DegDegM(:,3);
nedM = [northM,eastM,downM];
end %end of function Lla2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% function to difference latitude, longitude and altitude
% and provide an answer in NED coordinates in meters.
%
% Inputs: lla1DegDegM: mx3 matrix, [latitude(deg),longitude(deg),altitude(m)]
% lla2DegDegM: mx3 or 1x3, [latitude(deg),longitude(deg),altitude(m)]
% Output: nedM = lla1DegDegM - lla2DegDegM in NED coords (meters),
%
% Useful rules of thumb for quick conversions:
% 1e-5 (5th decimal place) of a degree of latitude approx= 1.1 meters
% 1e-5 of a degree of longitude approx= cos(latitude) * 1.1 meters
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if nargin<2, error('Two inputs arguments needed'),end
[m1,n1]=size(lla1DegDegM);
[m2,n2]=size(lla2DegDegM);
if m2==1
lla2DegDegM=ones(m1,1)*lla2DegDegM;
else
if m1~=m2,
error('Second input must have one row or same number of rows as first input')
end
end
if n1~=3 || n2~=3, error('Both input matrices must have 3 columns'),end
[xyz1M] = Lla2Xyz(lla1DegDegM);
[xyz2M] = Lla2Xyz(lla2DegDegM);
refXyz = (xyz1M+xyz2M)/2;
[llaDegDegM] = Xyz2Lla(refXyz);
northM = zeros(m1,1);
eastM = zeros(m1,1);
for i=1:m1
Ce2n = RotEcef2Ned(llaDegDegM(i,1),llaDegDegM(i,2));
v = Ce2n*(xyz1M(i,:)-xyz2M(i,:))';
northM(i)=v(1);
eastM(i)=v(2);
end
downM = -lla1DegDegM(:,3)+lla2DegDegM(:,3);
nedM = [northM,eastM,downM];
end %end of function Lla2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

108
opensource/Lla2Xyz.m → ProcessGNSSMeas/Lla2Xyz.m

@ -1,54 +1,54 @@
function [xyzM] = Lla2Xyz(llaDegDegM)
% [xyzM] = Lla2Xyz(llaDegDegM)
% Transform latitude, longitude, altitude to ECEF coordinates.
%
% input: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are in degrees, altitude in meters
% output: xyzM = [mx3] matrix of ECEF coordinates in meters
%
% See also Xyz2Lla
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(llaDegDegM,2)~=3
error('Input llaDegDegM must have three columns');
end
latDeg = llaDegDegM(:,1); lonDeg = llaDegDegM(:,2); altM = llaDegDegM(:,3);
%No rotation of longitude, by definition of ECEF
% Compute sines and cosines.
D2R = pi/180;
clat = cos(latDeg*D2R);
clon = cos(lonDeg*D2R);
slat = sin(latDeg*D2R);
slon = sin(lonDeg*D2R);
% Compute position vector in ECEF coordinates.
r0 = GpsConstants.EARTHSEMIMAJOR * ...
(sqrt(1.0 - GpsConstants.EARTHECCEN2 .* slat .* slat)).^(-1);
xM = (altM + r0) .* clat .* clon; % x coordinate
yM = (altM + r0) .* clat .* slon; % y coordinate
zM = (altM + r0 .* (1.0 - GpsConstants.EARTHECCEN2)).* slat;% z coordinate
[xyzM] = [xM,yM,zM];
end %end of function Lla2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [xyzM] = Lla2Xyz(llaDegDegM)
% [xyzM] = Lla2Xyz(llaDegDegM)
% Transform latitude, longitude, altitude to ECEF coordinates.
%
% input: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are in degrees, altitude in meters
% output: xyzM = [mx3] matrix of ECEF coordinates in meters
%
% See also Xyz2Lla
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(llaDegDegM,2)~=3
error('Input llaDegDegM must have three columns');
end
latDeg = llaDegDegM(:,1); lonDeg = llaDegDegM(:,2); altM = llaDegDegM(:,3);
%No rotation of longitude, by definition of ECEF
% Compute sines and cosines.
D2R = pi/180;
clat = cos(latDeg*D2R);
clon = cos(lonDeg*D2R);
slat = sin(latDeg*D2R);
slon = sin(lonDeg*D2R);
% Compute position vector in ECEF coordinates.
r0 = GpsConstants.EARTHSEMIMAJOR * ...
(sqrt(1.0 - GpsConstants.EARTHECCEN2 .* slat .* slat)).^(-1);
xM = (altM + r0) .* clat .* clon; % x coordinate
yM = (altM + r0) .* clat .* slon; % y coordinate
zM = (altM + r0 .* (1.0 - GpsConstants.EARTHECCEN2)).* slat;% z coordinate
[xyzM] = [xM,yM,zM];
end %end of function Lla2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/PlotAdr.m → ProcessGNSSMeas/PlotAdr.m

0
opensource/PlotAdrResids.m → ProcessGNSSMeas/PlotAdrResids.m

0
opensource/PlotCno.m → ProcessGNSSMeas/PlotCno.m

0
opensource/PlotPseudorangeRates.m → ProcessGNSSMeas/PlotPseudorangeRates.m

0
opensource/PlotPseudoranges.m → ProcessGNSSMeas/PlotPseudoranges.m

0
opensource/PlotPvt.m → ProcessGNSSMeas/PlotPvt.m

0
opensource/PlotPvtStates.m → ProcessGNSSMeas/PlotPvtStates.m

0
opensource/ProcessAdr.m → ProcessGNSSMeas/ProcessAdr.m

0
opensource/ProcessGnssMeas.m → ProcessGNSSMeas/ProcessGnssMeas.m

0
opensource/ProcessGnssMeasScript.m → ProcessGNSSMeas/ProcessGnssMeasScript.m

0
opensource/README.md → ProcessGNSSMeas/README.md

0
opensource/ReadGnssLogger.m → ProcessGNSSMeas/ReadGnssLogger.m

514
opensource/ReadRinexNav.m → ProcessGNSSMeas/ReadRinexNav.m

@ -1,257 +1,257 @@
function [gpsEph,iono] = ReadRinexNav(fileName)
% [gpsEph,iono] = ReadRinexNav(fileName)
%
% Read GPS ephemeris and iono data from an ASCII formatted RINEX 2.10 Nav file.
% Input:
% fileName - string containing name of RINEX formatted navigation data file
% Output:
% gpsEph: vector of ephemeris data, each element is an ephemeris structure
% structure order and orbit variable names follow RINEX 2.1 Table A4
% clock variable names af0, af1, af2 follow IS GPS 200
% gpsEph(i).PRN % SV PRN number
% gpsEph(i).Toc % Time of clock (seconds)
% gpsEph(i).af0 % SV clock bias (seconds)
% gpsEph(i).af1 % SV clock drift (sec/sec)
% gpsEph(i).af2 % SV clock drift rate (sec/sec2)
% gpsEph(i).IODE % Issue of data, ephemeris
% gpsEph(i).Crs % Sine harmonic correction to orbit radius (meters)
% gpsEph(i).Delta_n % Mean motion difference from computed value (radians/sec)
% gpsEph(i).M0 % Mean anomaly at reference time (radians)
% gpsEph(i).Cuc % Cosine harmonic correction to argument of lat (radians)
% gpsEph(i).e % Eccentricity (dimensionless)
% gpsEph(i).Cus % Sine harmonic correction to argument of latitude (radians)
% gpsEph(i).Asqrt % Square root of semi-major axis (meters^1/2)
% gpsEph(i).Toe % Reference time of ephemeris (seconds)
% gpsEph(i).Cic % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).OMEGA % Longitude of ascending node at weekly epoch (radians)
% gpsEph(i).Cis % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).i0 % Inclination angle at reference time (radians)
% gpsEph(i).Crc % Cosine harmonic correction to the orbit radius (meters)
% gpsEph(i).omega % Argument of perigee (radians)
% gpsEph(i).OMEGA_DOT% Rate of right ascension (radians/sec)
% gpsEph(i).IDOT % Rate of inclination angle (radians/sec)
% gpsEph(i).codeL2 % codes on L2 channel
% gpsEph(i).GPS_Week % GPS week (to go with Toe), (NOT Mod 1024)
% gpsEph(i).L2Pdata % L2 P data flag
% gpsEph(i).accuracy % SV user range accuracy (meters)
% gpsEph(i).health % Satellite health
% gpsEph(i).TGD % Group delay (seconds)
% gpsEph(i).IODC % Issue of Data, Clock
% gpsEph(i).ttx % Transmission time of message (seconds)
% gpsEph(i).Fit_interval %fit interval (hours), zero if not known
%
% iono: ionospheric parameter structure
% iono.alpha = [alpha0, alpha1, alpha2, alpha3]
% iono.beta = [ beta0, beta1, beta2, beta3]
% if iono data is not present in the Rinex file, iono is returned empty.
fidEph = fopen(fileName);
[numEph,numHdrLines] = countEph(fidEph);
%Now read from the begining again, looking for iono parameters
frewind(fidEph);
iono = readIono(fidEph,numHdrLines);
%initialize ephemeris structure array:
gpsEph = InitializeGpsEph;
gpsEph = repmat(gpsEph,1,numEph);
%now read each ephemeris into gpsEph(j)
%RINEX defines the format in terms of numbers of characters, so that's how we
%read it, e.g. "gpsEph(j).PRN = str2num(line(1:2));" and so on
for j = 1:numEph
line = fgetl(fidEph);
gpsEph(j).PRN = str2num(line(1:2));
%NOTE: we use str2num, not str2double, since str2num handles 'D' for exponent
%% get Toc (Rinex gives this as UTC time yy,mm,dd,hh,mm,ss)
year = str2num(line(3:6));
%convert year to a 4-digit year, this code is good to the year 2080.
%From 2080 RINEX 2.1 is ambiguous and shouldnt be used, because is has a
%2-digit year, and 100 years will have passed since the GPS Epoch.
if year < 80,
year = 2000+year;
else
year = 1900+year;
end
month = str2num(line(7:9));
day = str2num(line(10:12));
hour = str2num(line(13:15));
minute = str2num(line(16:18));
second = str2num(line(19:22));
%convert Toc to gpsTime
gpsTime = Utc2Gps([year,month,day,hour,minute,second]);
gpsEph(j).Toc = gpsTime(2);
%% get all other parameters
gpsEph(j).af0 = str2num(line(23:41));
gpsEph(j).af1 = str2num(line(42:60));
gpsEph(j).af2 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IODE = str2num(line(4:22));
gpsEph(j).Crs = str2num(line(23:41));
gpsEph(j).Delta_n = str2num(line(42:60));
gpsEph(j).M0 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).Cuc = str2num(line(4:22));
gpsEph(j).e = str2num(line(23:41));
gpsEph(j).Cus = str2num(line(42:60));
gpsEph(j).Asqrt = str2num(line(61:79));
line=fgetl(fidEph);
gpsEph(j).Toe = str2num(line(4:22));
gpsEph(j).Cic = str2num(line(23:41));
gpsEph(j).OMEGA = str2num(line(42:60));
gpsEph(j).Cis = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).i0 = str2num(line(4:22));
gpsEph(j).Crc = str2num(line(23:41));
gpsEph(j).omega = str2num(line(42:60));
gpsEph(j).OMEGA_DOT = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IDOT = str2num(line(4:22));
gpsEph(j).codeL2 = str2num(line(23:41));
gpsEph(j).GPS_Week = str2num(line(42:60));
gpsEph(j).L2Pdata = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).accuracy = str2num(line(4:22));
gpsEph(j).health = str2num(line(23:41));
gpsEph(j).TGD = str2num(line(42:60));
gpsEph(j).IODC = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).ttx = str2num(line(4:22));
gpsEph(j).Fit_interval = str2num(line(23:41));
end
fclose(fidEph);
end %end of function ReadRinexNav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [numEph,numHdrLines] = countEph(fidEph,fileName)
%utility function for ReadRinexNav
%Read past the header, and then read to the end, counting ephemerides:
numHdrLines = 0;
bFoundHeader = false;
while ~bFoundHeader %Read past the header
numHdrLines = numHdrLines+1;
line = fgetl(fidEph);
if ~ischar(line), break, end
k = strfind(line,'END OF HEADER');
if ~isempty(k),
bFoundHeader = true;
break
end
end
if ~bFoundHeader
error('Error reading file: %s\nExpected RINEX header not found',fileName);
end
%Now read to the end of the file
numEph = -1;
while 1
numEph = numEph+1;
line = fgetl(fidEph);
if line == -1,
break;
elseif length(line)~=79
%use this opportunity to check line is the right length
%because in the rest of ReadRinexNav we depend on line being this length
error('Incorrect line length encountered in RINEX file');
end
end;
%check that we read the expected number of lines:
if rem(numEph,8)~=0
error('Number of nav lines in %s should be divisible by 8',fileName);
end
numEph = numEph/8; %8 lines per ephemeris
end %end of function countEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function iono = readIono(fidEph,numHdrLines)
%utility function to read thru the header lines, and find iono parameters
iono = []; %return empty if iono not found
bIonoAlpha=false; bIonoBeta=false;
for i = 1:numHdrLines,
line = fgetl(fidEph);
% Look for iono parameters, and read them in
if ~isempty(strfind(line,'ION ALPHA')) %line contains iono alpha parameters
ii = strfind(line,'ION ALPHA');
iono.alpha=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono alpha
bIonoAlpha = (length(iono.alpha)==4);
end
if ~isempty(strfind(line,'ION BETA'))%line contains the iono beta parameters
ii = strfind(line,'ION BETA');
iono.beta=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono beta
bIonoBeta = (length(iono.beta)==4);
end
end
if ~(bIonoAlpha && bIonoBeta)
%if we didn't get both alpha and beta iono correctly, then return empty iono
iono=[];
end
end %end of function readIono
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gpsEph = InitializeGpsEph
%utility function to initialize the ephemeris structure
gpsEph.PRN = 0;
gpsEph.Toc = 0;
gpsEph.af0 = 0;
gpsEph.af1 = 0;
gpsEph.af2 = 0;
gpsEph.IODE = 0;
gpsEph.Crs = 0;
gpsEph.Delta_n = 0;
gpsEph.M0 = 0;
gpsEph.Cuc = 0;
gpsEph.e = 0;
gpsEph.Cus = 0;
gpsEph.Asqrt = 0;
gpsEph.Toe = 0;
gpsEph.Cic = 0;
gpsEph.OMEGA = 0;
gpsEph.Cis = 0;
gpsEph.i0 = 0;
gpsEph.Crc = 0;
gpsEph.omega = 0;
gpsEph.OMEGA_DOT = 0;
gpsEph.IDOT = 0;
gpsEph.codeL2 = 0;
gpsEph.GPS_Week = 0;
gpsEph.L2Pdata = 0;
gpsEph.accuracy = 0;
gpsEph.health = 0;
gpsEph.TGD = 0;
gpsEph.IODC = 0;
gpsEph.ttx = 0;
gpsEph.Fit_interval= 0;
end %end of function InitializeEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [gpsEph,iono] = ReadRinexNav(fileName)
% [gpsEph,iono] = ReadRinexNav(fileName)
%
% Read GPS ephemeris and iono data from an ASCII formatted RINEX 2.10 Nav file.
% Input:
% fileName - string containing name of RINEX formatted navigation data file
% Output:
% gpsEph: vector of ephemeris data, each element is an ephemeris structure
% structure order and orbit variable names follow RINEX 2.1 Table A4
% clock variable names af0, af1, af2 follow IS GPS 200
% gpsEph(i).PRN % SV PRN number
% gpsEph(i).Toc % Time of clock (seconds)
% gpsEph(i).af0 % SV clock bias (seconds)
% gpsEph(i).af1 % SV clock drift (sec/sec)
% gpsEph(i).af2 % SV clock drift rate (sec/sec2)
% gpsEph(i).IODE % Issue of data, ephemeris
% gpsEph(i).Crs % Sine harmonic correction to orbit radius (meters)
% gpsEph(i).Delta_n % Mean motion difference from computed value (radians/sec)
% gpsEph(i).M0 % Mean anomaly at reference time (radians)
% gpsEph(i).Cuc % Cosine harmonic correction to argument of lat (radians)
% gpsEph(i).e % Eccentricity (dimensionless)
% gpsEph(i).Cus % Sine harmonic correction to argument of latitude (radians)
% gpsEph(i).Asqrt % Square root of semi-major axis (meters^1/2)
% gpsEph(i).Toe % Reference time of ephemeris (seconds)
% gpsEph(i).Cic % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).OMEGA % Longitude of ascending node at weekly epoch (radians)
% gpsEph(i).Cis % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).i0 % Inclination angle at reference time (radians)
% gpsEph(i).Crc % Cosine harmonic correction to the orbit radius (meters)
% gpsEph(i).omega % Argument of perigee (radians)
% gpsEph(i).OMEGA_DOT% Rate of right ascension (radians/sec)
% gpsEph(i).IDOT % Rate of inclination angle (radians/sec)
% gpsEph(i).codeL2 % codes on L2 channel
% gpsEph(i).GPS_Week % GPS week (to go with Toe), (NOT Mod 1024)
% gpsEph(i).L2Pdata % L2 P data flag
% gpsEph(i).accuracy % SV user range accuracy (meters)
% gpsEph(i).health % Satellite health
% gpsEph(i).TGD % Group delay (seconds)
% gpsEph(i).IODC % Issue of Data, Clock
% gpsEph(i).ttx % Transmission time of message (seconds)
% gpsEph(i).Fit_interval %fit interval (hours), zero if not known
%
% iono: ionospheric parameter structure
% iono.alpha = [alpha0, alpha1, alpha2, alpha3]
% iono.beta = [ beta0, beta1, beta2, beta3]
% if iono data is not present in the Rinex file, iono is returned empty.
fidEph = fopen(fileName);
[numEph,numHdrLines] = countEph(fidEph);
%Now read from the begining again, looking for iono parameters
frewind(fidEph);
iono = readIono(fidEph,numHdrLines);
%initialize ephemeris structure array:
gpsEph = InitializeGpsEph;
gpsEph = repmat(gpsEph,1,numEph);
%now read each ephemeris into gpsEph(j)
%RINEX defines the format in terms of numbers of characters, so that's how we
%read it, e.g. "gpsEph(j).PRN = str2num(line(1:2));" and so on
for j = 1:numEph
line = fgetl(fidEph);
gpsEph(j).PRN = str2num(line(1:2));
%NOTE: we use str2num, not str2double, since str2num handles 'D' for exponent
%% get Toc (Rinex gives this as UTC time yy,mm,dd,hh,mm,ss)
year = str2num(line(3:6));
%convert year to a 4-digit year, this code is good to the year 2080.
%From 2080 RINEX 2.1 is ambiguous and shouldnt be used, because is has a
%2-digit year, and 100 years will have passed since the GPS Epoch.
if year < 80,
year = 2000+year;
else
year = 1900+year;
end
month = str2num(line(7:9));
day = str2num(line(10:12));
hour = str2num(line(13:15));
minute = str2num(line(16:18));
second = str2num(line(19:22));
%convert Toc to gpsTime
gpsTime = Utc2Gps([year,month,day,hour,minute,second]);
gpsEph(j).Toc = gpsTime(2);
%% get all other parameters
gpsEph(j).af0 = str2num(line(23:41));
gpsEph(j).af1 = str2num(line(42:60));
gpsEph(j).af2 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IODE = str2num(line(4:22));
gpsEph(j).Crs = str2num(line(23:41));
gpsEph(j).Delta_n = str2num(line(42:60));
gpsEph(j).M0 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).Cuc = str2num(line(4:22));
gpsEph(j).e = str2num(line(23:41));
gpsEph(j).Cus = str2num(line(42:60));
gpsEph(j).Asqrt = str2num(line(61:79));
line=fgetl(fidEph);
gpsEph(j).Toe = str2num(line(4:22));
gpsEph(j).Cic = str2num(line(23:41));
gpsEph(j).OMEGA = str2num(line(42:60));
gpsEph(j).Cis = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).i0 = str2num(line(4:22));
gpsEph(j).Crc = str2num(line(23:41));
gpsEph(j).omega = str2num(line(42:60));
gpsEph(j).OMEGA_DOT = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IDOT = str2num(line(4:22));
gpsEph(j).codeL2 = str2num(line(23:41));
gpsEph(j).GPS_Week = str2num(line(42:60));
gpsEph(j).L2Pdata = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).accuracy = str2num(line(4:22));
gpsEph(j).health = str2num(line(23:41));
gpsEph(j).TGD = str2num(line(42:60));
gpsEph(j).IODC = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).ttx = str2num(line(4:22));
gpsEph(j).Fit_interval = str2num(line(23:41));
end
fclose(fidEph);
end %end of function ReadRinexNav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [numEph,numHdrLines] = countEph(fidEph,fileName)
%utility function for ReadRinexNav
%Read past the header, and then read to the end, counting ephemerides:
numHdrLines = 0;
bFoundHeader = false;
while ~bFoundHeader %Read past the header
numHdrLines = numHdrLines+1;
line = fgetl(fidEph);
if ~ischar(line), break, end
k = strfind(line,'END OF HEADER');
if ~isempty(k),
bFoundHeader = true;
break
end
end
if ~bFoundHeader
error('Error reading file: %s\nExpected RINEX header not found',fileName);
end
%Now read to the end of the file
numEph = -1;
while 1
numEph = numEph+1;
line = fgetl(fidEph);
if line == -1,
break;
elseif length(line)~=79
%use this opportunity to check line is the right length
%because in the rest of ReadRinexNav we depend on line being this length
error('Incorrect line length encountered in RINEX file');
end
end;
%check that we read the expected number of lines:
if rem(numEph,8)~=0
error('Number of nav lines in %s should be divisible by 8',fileName);
end
numEph = numEph/8; %8 lines per ephemeris
end %end of function countEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function iono = readIono(fidEph,numHdrLines)
%utility function to read thru the header lines, and find iono parameters
iono = []; %return empty if iono not found
bIonoAlpha=false; bIonoBeta=false;
for i = 1:numHdrLines,
line = fgetl(fidEph);
% Look for iono parameters, and read them in
if ~isempty(strfind(line,'ION ALPHA')) %line contains iono alpha parameters
ii = strfind(line,'ION ALPHA');
iono.alpha=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono alpha
bIonoAlpha = (length(iono.alpha)==4);
end
if ~isempty(strfind(line,'ION BETA'))%line contains the iono beta parameters
ii = strfind(line,'ION BETA');
iono.beta=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono beta
bIonoBeta = (length(iono.beta)==4);
end
end
if ~(bIonoAlpha && bIonoBeta)
%if we didn't get both alpha and beta iono correctly, then return empty iono
iono=[];
end
end %end of function readIono
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gpsEph = InitializeGpsEph
%utility function to initialize the ephemeris structure
gpsEph.PRN = 0;
gpsEph.Toc = 0;
gpsEph.af0 = 0;
gpsEph.af1 = 0;
gpsEph.af2 = 0;
gpsEph.IODE = 0;
gpsEph.Crs = 0;
gpsEph.Delta_n = 0;
gpsEph.M0 = 0;
gpsEph.Cuc = 0;
gpsEph.e = 0;
gpsEph.Cus = 0;
gpsEph.Asqrt = 0;
gpsEph.Toe = 0;
gpsEph.Cic = 0;
gpsEph.OMEGA = 0;
gpsEph.Cis = 0;
gpsEph.i0 = 0;
gpsEph.Crc = 0;
gpsEph.omega = 0;
gpsEph.OMEGA_DOT = 0;
gpsEph.IDOT = 0;
gpsEph.codeL2 = 0;
gpsEph.GPS_Week = 0;
gpsEph.L2Pdata = 0;
gpsEph.accuracy = 0;
gpsEph.health = 0;
gpsEph.TGD = 0;
gpsEph.IODC = 0;
gpsEph.ttx = 0;
gpsEph.Fit_interval= 0;
end %end of function InitializeEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

120
opensource/RotEcef2Ned.m → ProcessGNSSMeas/RotEcef2Ned.m

@ -1,60 +1,60 @@
function Re2n = RotEcef2Ned(latDeg, lonDeg)
% Re2n = RotEcef2Ned(latDeg, lonDeg)
% Rotation matrix to convert an ECEF vector to
% North, East, Down coordinates, and vice-versa
%
% inputs: latDeg, lonDeg (degrees)
% output: Re2n, 3x3 unitary rotation matrix =
% [-sin(lat)*cos(lon), -sin(lat)*sin(lon), cos(lat);
% -sin(lon), cos(lon), 0 ;
% -cos(lat)*cos(lon), -cos(lat)*sin(lon),-sin(lat)]
%
% Example: vNed = Re2n*vEcef,
% Re2n'*vNed = vEcef
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%CHECK INPUTS
if any(size(latDeg)~=[1,1]) || any(size(lonDeg)~=[1,1])
error('Inputs latDeg, lonDeg must be scalars')
end
D2R = pi/180; %degrees to radians scale factor
latRad=D2R*latDeg(:); lonRad=D2R*lonDeg(:);
clat = cos(latRad);
slat = sin(latRad);
clon = cos(lonRad);
slon = sin(lonRad);
Re2n = zeros(3,3);
Re2n(1,1) = -slat.*clon;
Re2n(1,2) = -slat.*slon;
Re2n(1,3) = clat;
Re2n(2,1) = -slon;
Re2n(2,2) = clon;
Re2n(2,3) = 0;
Re2n(3,1) = -clat.*clon;
Re2n(3,2) = -clat.*slon;
Re2n(3,3) = -slat;
end %end of function RotEcef2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function Re2n = RotEcef2Ned(latDeg, lonDeg)
% Re2n = RotEcef2Ned(latDeg, lonDeg)
% Rotation matrix to convert an ECEF vector to
% North, East, Down coordinates, and vice-versa
%
% inputs: latDeg, lonDeg (degrees)
% output: Re2n, 3x3 unitary rotation matrix =
% [-sin(lat)*cos(lon), -sin(lat)*sin(lon), cos(lat);
% -sin(lon), cos(lon), 0 ;
% -cos(lat)*cos(lon), -cos(lat)*sin(lon),-sin(lat)]
%
% Example: vNed = Re2n*vEcef,
% Re2n'*vNed = vEcef
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%CHECK INPUTS
if any(size(latDeg)~=[1,1]) || any(size(lonDeg)~=[1,1])
error('Inputs latDeg, lonDeg must be scalars')
end
D2R = pi/180; %degrees to radians scale factor
latRad=D2R*latDeg(:); lonRad=D2R*lonDeg(:);
clat = cos(latRad);
slat = sin(latRad);
clon = cos(lonRad);
slon = sin(lonRad);
Re2n = zeros(3,3);
Re2n(1,1) = -slat.*clon;
Re2n(1,2) = -slat.*slon;
Re2n(1,3) = clat;
Re2n(2,1) = -slon;
Re2n(2,2) = clon;
Re2n(2,3) = 0;
Re2n(3,1) = -clat.*clon;
Re2n(3,2) = -clat.*slon;
Re2n(3,3) = -slat;
end %end of function RotEcef2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/SetDataFilter.m → ProcessGNSSMeas/SetDataFilter.m

222
opensource/Utc2Gps.m → ProcessGNSSMeas/Utc2Gps.m

@ -1,111 +1,111 @@
function [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% Convert the UTC date and time to GPS week & seconds
%
% Inputs:
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Outputs:
% gpsTime: [mx2] matrix [gpsWeek, gpsSeconds],
% gpsWeek = number of weeks since GPS epoch
% gpsSeconds = number of seconds into gpsWeek,
% fctSeconds: full cycle time = seconds since GPS epoch (1980/01/06 00:00 UTC)
% Other functions needed: JulianDay.m, LeapSeconds.m
%initialize outputs
gpsTime=[];
fctSeconds=[];
[bOk]=CheckUtcTimeInputs(utcTime);
if ~bOk
return
end
HOURSEC = 3600; MINSEC = 60;
daysSinceEpoch = floor(JulianDay(utcTime) - GpsConstants.GPSEPOCHJD);
gpsWeek = fix(daysSinceEpoch/7);
dayofweek = rem(daysSinceEpoch,7);
% calculate the number of seconds since Sunday at midnight:
gpsSeconds = dayofweek*GpsConstants.DAYSEC + utcTime(:,4)*HOURSEC + ...
utcTime(:,5)*MINSEC + utcTime(:,6);
gpsWeek = gpsWeek + fix(gpsSeconds/GpsConstants.WEEKSEC);
gpsSeconds = rem(gpsSeconds,GpsConstants.WEEKSEC);
% now add leapseconds
leapSecs = LeapSeconds(utcTime);
fctSeconds = gpsWeek(:)*GpsConstants.WEEKSEC + gpsSeconds(:) + leapSecs(:);
% when a leap second happens, utc time stands still for one second,
% so gps seconds get further ahead, so we add leapsecs in going to gps time
gpsWeek = fix(fctSeconds/GpsConstants.WEEKSEC);
iZ = gpsWeek==0;
gpsSeconds(iZ) = fctSeconds(iZ); %set gpsSeconds directly, because rem(x,0)=NaN
gpsSeconds(~iZ) = rem(fctSeconds(~iZ),gpsWeek(~iZ)*GpsConstants.WEEKSEC);
gpsTime=[gpsWeek,gpsSeconds];
assert(all(fctSeconds==gpsWeek*GpsConstants.WEEKSEC+gpsSeconds),...
'Error in computing gpsWeek, gpsSeconds');
end %end of function Utc2Gps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bOk] = CheckUtcTimeInputs(utcTime)
%utility function for Utc2Gps
%check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
%check that year, month & day are integers
x = utcTime(:,1:3);
if any(any( (x-fix(x)) ~= 0 ))
error('year,month & day must be integers')
end
%check that year is in valid range
if ( any(utcTime(:,1)<1980) || any(utcTime(:,1)>2099) )
error('year must have values in the range: [1980:2099]')
end
%check validity of month, day and time
if (any(utcTime(:,2))<1 || any(utcTime(:,2))>12 )
error('The month in utcTime must be a number in the set [1:12]')
end
if (any(utcTime(:,3)<1) || any(utcTime(:,3)>31))
error('The day in utcTime must be a number in the set [1:31]')
end
if (any(utcTime(:,4)<0) || any(utcTime(:,4)>=24))
error('The hour in utcTime must be in the range [0,24)')
end
if (any(utcTime(:,5)<0) || any(utcTime(:,5)>=60))
error('The minutes in utcTime must be in the range [0,60)')
end
if (any(utcTime(:,6)<0) || any(utcTime(:,6)>60))
%Note: seconds can equal 60 exactly, on the second a leap second is added
error('The seconds in utcTime must be in the range [0,60]')
end
bOk = true;
end %end of function CheckUtcTimeInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% Convert the UTC date and time to GPS week & seconds
%
% Inputs:
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Outputs:
% gpsTime: [mx2] matrix [gpsWeek, gpsSeconds],
% gpsWeek = number of weeks since GPS epoch
% gpsSeconds = number of seconds into gpsWeek,
% fctSeconds: full cycle time = seconds since GPS epoch (1980/01/06 00:00 UTC)
% Other functions needed: JulianDay.m, LeapSeconds.m
%initialize outputs
gpsTime=[];
fctSeconds=[];
[bOk]=CheckUtcTimeInputs(utcTime);
if ~bOk
return
end
HOURSEC = 3600; MINSEC = 60;
daysSinceEpoch = floor(JulianDay(utcTime) - GpsConstants.GPSEPOCHJD);
gpsWeek = fix(daysSinceEpoch/7);
dayofweek = rem(daysSinceEpoch,7);
% calculate the number of seconds since Sunday at midnight:
gpsSeconds = dayofweek*GpsConstants.DAYSEC + utcTime(:,4)*HOURSEC + ...
utcTime(:,5)*MINSEC + utcTime(:,6);
gpsWeek = gpsWeek + fix(gpsSeconds/GpsConstants.WEEKSEC);
gpsSeconds = rem(gpsSeconds,GpsConstants.WEEKSEC);
% now add leapseconds
leapSecs = LeapSeconds(utcTime);
fctSeconds = gpsWeek(:)*GpsConstants.WEEKSEC + gpsSeconds(:) + leapSecs(:);
% when a leap second happens, utc time stands still for one second,
% so gps seconds get further ahead, so we add leapsecs in going to gps time
gpsWeek = fix(fctSeconds/GpsConstants.WEEKSEC);
iZ = gpsWeek==0;
gpsSeconds(iZ) = fctSeconds(iZ); %set gpsSeconds directly, because rem(x,0)=NaN
gpsSeconds(~iZ) = rem(fctSeconds(~iZ),gpsWeek(~iZ)*GpsConstants.WEEKSEC);
gpsTime=[gpsWeek,gpsSeconds];
assert(all(fctSeconds==gpsWeek*GpsConstants.WEEKSEC+gpsSeconds),...
'Error in computing gpsWeek, gpsSeconds');
end %end of function Utc2Gps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bOk] = CheckUtcTimeInputs(utcTime)
%utility function for Utc2Gps
%check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
%check that year, month & day are integers
x = utcTime(:,1:3);
if any(any( (x-fix(x)) ~= 0 ))
error('year,month & day must be integers')
end
%check that year is in valid range
if ( any(utcTime(:,1)<1980) || any(utcTime(:,1)>2099) )
error('year must have values in the range: [1980:2099]')
end
%check validity of month, day and time
if (any(utcTime(:,2))<1 || any(utcTime(:,2))>12 )
error('The month in utcTime must be a number in the set [1:12]')
end
if (any(utcTime(:,3)<1) || any(utcTime(:,3)>31))
error('The day in utcTime must be a number in the set [1:31]')
end
if (any(utcTime(:,4)<0) || any(utcTime(:,4)>=24))
error('The hour in utcTime must be in the range [0,24)')
end
if (any(utcTime(:,5)<0) || any(utcTime(:,5)>=60))
error('The minutes in utcTime must be in the range [0,60)')
end
if (any(utcTime(:,6)<0) || any(utcTime(:,6)>60))
%Note: seconds can equal 60 exactly, on the second a leap second is added
error('The seconds in utcTime must be in the range [0,60]')
end
bOk = true;
end %end of function CheckUtcTimeInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

0
opensource/WlsPvt.m → ProcessGNSSMeas/WlsPvt.m

164
opensource/Xyz2Lla.m → ProcessGNSSMeas/Xyz2Lla.m

@ -1,82 +1,82 @@
function [llaDegDegM] = Xyz2Lla(xyzM)
% [llaDegDegM] = Xyz2Lla(xyzM)
%
% Transform Earth-Centered-Earth-Fixed x,y,z, coordinates to lat,lon,alt.
% Input: xyzM = [mx3] matrix of ECEF coordinates in meters
% Output: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are returned in degrees and altitude in meters
%
% See also Lla2Xyz
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(xyzM,2)~=3
error('Input xyzM must have three columns');
end
% algorithm: Hoffman-Wellenhof, Lichtenegger & Collins "GPS Theory & Practice"
R2D = 180/pi;
%if x and y ecef positions are both zero then lla is undefined
iZero = ( xyzM(:,1)==0 & xyzM(:,2)==0);
xyzM(iZero,:) = NaN; %set to NaN, so lla will also be NaN
xM = xyzM(:,1); yM = xyzM(:,2); zM = xyzM(:,3);
%following algorithm from Hoffman-Wellenhof, et al. "GPS Theory & Practice":
a = GpsConstants.EARTHSEMIMAJOR;
a2 = a^2;
b2 = a2*(1-GpsConstants.EARTHECCEN2);
b = sqrt(b2);
e2 = GpsConstants.EARTHECCEN2;
ep2 = (a2-b2)/b2;
p=sqrt(xM.^2 + yM.^2);
% two sides and hypotenuse of right angle triangle with one angle = theta:
s1 = zM*a;
s2 = p*b;
h = sqrt(s1.^2 + s2.^2);
sin_theta = s1./h;
cos_theta = s2./h;
%theta = atan(s1./s2);
% two sides and hypotenuse of right angle triangle with one angle = lat:
s1 = zM+ep2*b*(sin_theta.^3);
s2 = p-a*e2.*(cos_theta.^3);
h = sqrt(s1.^2 + s2.^2);
tan_lat = s1./s2;
sin_lat = s1./h;
cos_lat = s2./h;
latDeg = atan(tan_lat);
latDeg = latDeg*R2D;
N = a2*((a2*(cos_lat.^2) + b2*(sin_lat.^2)).^(-0.5));
altM = p./cos_lat - N;
% rotate longitude to where it would be for a fixed point in ECI
%lonDeg = atan2(yM,xM) - GpsConstants.WE*deltaTime;
lonDeg = atan2(yM,xM); %since deltaTime = 0 for ECEF
lonDeg = rem(lonDeg,2*pi)*R2D;
if (lonDeg>180)
lonDeg = lonDeg-360;
end
llaDegDegM = [latDeg, lonDeg,altM];
end %end of function Xyz2Lla
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function [llaDegDegM] = Xyz2Lla(xyzM)
% [llaDegDegM] = Xyz2Lla(xyzM)
%
% Transform Earth-Centered-Earth-Fixed x,y,z, coordinates to lat,lon,alt.
% Input: xyzM = [mx3] matrix of ECEF coordinates in meters
% Output: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are returned in degrees and altitude in meters
%
% See also Lla2Xyz
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(xyzM,2)~=3
error('Input xyzM must have three columns');
end
% algorithm: Hoffman-Wellenhof, Lichtenegger & Collins "GPS Theory & Practice"
R2D = 180/pi;
%if x and y ecef positions are both zero then lla is undefined
iZero = ( xyzM(:,1)==0 & xyzM(:,2)==0);
xyzM(iZero,:) = NaN; %set to NaN, so lla will also be NaN
xM = xyzM(:,1); yM = xyzM(:,2); zM = xyzM(:,3);
%following algorithm from Hoffman-Wellenhof, et al. "GPS Theory & Practice":
a = GpsConstants.EARTHSEMIMAJOR;
a2 = a^2;
b2 = a2*(1-GpsConstants.EARTHECCEN2);
b = sqrt(b2);
e2 = GpsConstants.EARTHECCEN2;
ep2 = (a2-b2)/b2;
p=sqrt(xM.^2 + yM.^2);
% two sides and hypotenuse of right angle triangle with one angle = theta:
s1 = zM*a;
s2 = p*b;
h = sqrt(s1.^2 + s2.^2);
sin_theta = s1./h;
cos_theta = s2./h;
%theta = atan(s1./s2);
% two sides and hypotenuse of right angle triangle with one angle = lat:
s1 = zM+ep2*b*(sin_theta.^3);
s2 = p-a*e2.*(cos_theta.^3);
h = sqrt(s1.^2 + s2.^2);
tan_lat = s1./s2;
sin_lat = s1./h;
cos_lat = s2./h;
latDeg = atan(tan_lat);
latDeg = latDeg*R2D;
N = a2*((a2*(cos_lat.^2) + b2*(sin_lat.^2)).^(-0.5));
altM = p./cos_lat - N;
% rotate longitude to where it would be for a fixed point in ECI
%lonDeg = atan2(yM,xM) - GpsConstants.WE*deltaTime;
lonDeg = atan2(yM,xM); %since deltaTime = 0 for ECEF
lonDeg = rem(lonDeg,2*pi)*R2D;
if (lonDeg>180)
lonDeg = lonDeg-360;
end
llaDegDegM = [latDeg, lonDeg,altM];
end %end of function Xyz2Lla
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

3352
ProcessGNSSMeas/demoFiles/hour1820.16n

File diff suppressed because it is too large Load Diff

0
opensource/demoFiles/hour2570.16n → ProcessGNSSMeas/demoFiles/hour2570.16n

3478
ProcessGNSSMeas/demoFiles/raw.csv

File diff suppressed because it is too large Load Diff

0
opensource/demoFiles/workshop_trials01.txt → ProcessGNSSMeas/demoFiles/workshop_trials01.txt

Loading…
Cancel
Save