SPP

.vscode/tasks.json

@@ -0,0 +1,28 @@
+{
+    "tasks": [
+        {
+            "type": "cppbuild",
+            "label": "C/C++: gcc.exe 生成活动文件",
+            "command": "C:\\Users\\34716\\mingw64\\bin\\gcc.exe",
+            "args": [
+                "-fdiagnostics-color=always",
+                "-g",
+                "${file}",
+                "-o",
+                "${fileDirname}\\${fileBasenameNoExtension}.exe"
+            ],
+            "options": {
+                "cwd": "${fileDirname}"
+            },
+            "problemMatcher": [
+                "$gcc"
+            ],
+            "group": {
+                "kind": "build",
+                "isDefault": true
+            },
+            "detail": "调试器生成的任务。"
+        }
+    ],
+    "version": "2.0.0"
+}

2gpst.h

@@ -0,0 +1,26 @@
+#ifndef _2GPST_H_
+#define _2GPST_H_
+
+#include"read.h"
+
+double TimetoGPSsec(int y, int m, int d, double h, int min, double sec)
+{
+	if (m > 2)
+	{
+		y = y;
+		m = m;
+	}
+	if (m <= 2)
+	{
+		y = y - 1;
+		m = m + 12;
+	}	
+	h = h + min/60.0 + sec/3600.0;
+	double JD = (int)(365.25 * y) + (int)(30.6001 * (m + 1)) + d + h/24.0 + 1720981.5;
+	double MJD = JD - 2400000.5; 
+	int gpsweek = (int)((MJD - 44244) / 7);
+	double secofweek = (MJD - 44244.0 - gpsweek * 7.0) * 86400.0;
+	return secofweek;
+}
+
+#endif

CppProperties.json

@@ -0,0 +1,21 @@
+{
+  "configurations": [
+    {
+      "inheritEnvironments": [
+        "msvc_x64"
+      ],
+      "name": "x64-Debug",
+      "includePath": [
+        "${env.INCLUDE}",
+        "${workspaceRoot}\\**"
+      ],
+      "defines": [
+        "WIN32",
+        "_DEBUG",
+        "UNICODE",
+        "_UNICODE"
+      ],
+      "intelliSenseMode": "windows-msvc-x64"
+    }
+  ]
+}

Draw the Track[by PyGMT]/.idea/.gitignore

@@ -0,0 +1,3 @@
+# 默认忽略的文件
+/shelf/
+/workspace.xml

Draw the Track[by PyGMT]/.idea/GMT4SAT.iml

@@ -0,0 +1,15 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="jdk" jdkName="Python 3.8 (pygmt)" jdkType="Python SDK" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+  <component name="PyDocumentationSettings">
+    <option name="format" value="PLAIN" />
+    <option name="myDocStringFormat" value="Plain" />
+  </component>
+  <component name="TestRunnerService">
+    <option name="PROJECT_TEST_RUNNER" value="py.test" />
+  </component>
+</module>

Draw the Track[by PyGMT]/.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

Draw the Track[by PyGMT]/.idea/misc.xml

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="Black">
+    <option name="sdkName" value="Python 3.12" />
+  </component>
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (pygmt)" project-jdk-type="Python SDK" />
+</project>

Draw the Track[by PyGMT]/.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/GMT4SAT.iml" filepath="$PROJECT_DIR$/.idea/GMT4SAT.iml" />
+    </modules>
+  </component>
+</project>

Draw the Track[by PyGMT]/pyGMT.py

@@ -0,0 +1,98 @@
+import os
+
+from pygmt import project
+
+os.environ["GMT_LIBRARY_PATH"] = "C:/programs/gmt6/bin"
+
+import pygmt
+
+fig = pygmt.Figure()
+fig.coast(region="d",projection="Cyl_stere/12c",land="white",water="#C6E2FF",borders="1/0.2p",shorelines="1/0.2p",frame="ag")
+
+f = open("LLA_result_for_PyGMT.txt","r")
+f.readline()
+string = f.readline()
+epoch = 0
+while epoch < 2880:
+    flag = string[0:1]
+    epoch = int(string[1:5])
+    if flag == ">":
+        gpsnum = 0
+        while gpsnum != 36:
+            string = f.readline()
+            flag = string[0:1]
+            if flag =="G":
+                gpsnum = gpsnum + 1
+                sPRN = int(string[1:3])
+                B = float(string[10:20])
+                L = float(string[21:31])
+                if sPRN == 1:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 2:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 3:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 4:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 5:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 6:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 7:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 8:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 9:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 10:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 11:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 12:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 13:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 14:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 15:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 16:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 17:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 18:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 19:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 20:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 21:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 22:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 23:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 24:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 25:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 26:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 27:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 28:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 29:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 30:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 31:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+                if sPRN == 32:
+                    fig.plot(x = L,y = B,style = "c0.015c",fill = "blue",pen = "blue")
+            if flag !="G":
+                break
+    elif flag != ">":
+        break
+
+fig.show()

Observation_Station_Site-Solving_by_Matlab/Station_Site_Solving.asv

@@ -0,0 +1,170 @@
+% 打开坐标文件
+file=fopen('LLA_result_for_Matlab.txt');       
+sPRN = [];
+satX = [];
+satY = [];
+satZ = [];
+P_t1 = [];
+P_t2 = [];
+H    = [];
+TGD  = [];
+Dt   = [];
+sat_num = [];
+
+C_V    = 299792458;
+f1     =   1575.42;
+f2     =   1227.60;
+satnum =         0;
+eponum =         0;
+
+% 设置测站粗略坐标和精确的经纬高数据
+apX = -4052052.8310;
+apY =  4212835.9820;
+apZ = -2545104.4760;
+
+apB = -0.4131213;%纬度
+apL =  2.3367431;%经度
+apH =    603.235;
+
+% 逐行读取文件
+while ~feof(file)
+        line = fgets(file);
+        flag_epo = extractBetween(line,1 ,1 );
+%   判断历元标识
+    if strcmp(flag_epo, ">") == 1
+         eponum = str2double(extractBetween(line,2 ,5 ));
+       while eponum <= 2880 
+               flag_sat = extractBetween(line,1 ,1);
+%          判断卫星标识，读取卫星PRN号和XYZ坐标
+           if strcmp(flag_sat, "G") == 1
+                satnum = satnum + 1;
+                P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16               ,25 + 16 + 16 + 16 ));
+                if P_t1(eponum, satnum) ~= 0
+                    sPRN(eponum, satnum) = str2double(extractBetween(line,2 ,3 ));
+                    satX(eponum, satnum) = str2double(extractBetween(line,11                                   ,25 ));
+                    satY(eponum, satnum) = str2double(extractBetween(line,11 + 16                              ,25 + 16 ));
+                    satZ(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16                         ,25 + 16 + 16 ));
+                    P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16                    ,25 + 16 + 16 + 16 ));
+                    P_t2(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16               ,25 + 16 + 16 + 16 + 16));
+                    H(eponum, satnum)    = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16          ,25 + 16 + 16 + 16 + 16 + 16));
+                    TGD(eponum, satnum)  = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16     ,25 + 16 + 16 + 16 + 16 + 16 + 16));
+                    Dt(eponum, satnum)   = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16 + 16,25 + 16 + 16 + 16 + 16 + 16 + 16 + 17));
+                    line = fgets(file);
+                end
+                if P_t1(eponum, satnum) == 0
+                    satnum = satnum - 1;
+                    line = fgets(file);
+                end
+           end
+%          开头字符不是G，但是是>,判断为一个历元的起点，进行数据读取
+           if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 1
+                satnum = 0;
+                line = fgets(file);
+                continue;
+           end
+%          开头字符既不是G，也不是>,判断为历元结束标志，跳出循环
+           if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 0 && ~feof(file)
+                sat_num(eponum) = satnum;
+                satnum = 0;
+                break;
+           end
+       end
+    end
+    if strcmp(flag_epo, ">") == 0
+        continue;
+    end
+end
+fclose(file);
+
+% 平差计算测站位置和接收机钟差解算
+for i = 1:1:eponum
+    if sat_num(i) >= 4 
+        X = zeros(4, 1);
+        P = zeros(sat_num(i), sat_num(i));
+        B = zeros(sat_num(i), 4);
+        l = zeros(sat_num(i), 1);
+        X_R = zeros(3, 1);
+        
+        DeltaXYZ = zeros(3, 1);
+%       每个历元构成系数矩阵
+        for j = 1:1:sat_num(i)
+            R = sqrt((satX(i, j)-apX)^2+(satY(i, j)-apY)^2+(satZ(i, j)-apZ)^2);
+            B(j, 1) = -((satX(i, j) - apX) / R);
+            B(j, 2) = -((satY(i, j) - apY) / R);
+            B(j, 3) = -((satZ(i, j) - apZ) / R);
+            B(j, 4) = 1;
+
+            dtrop = 2.47/sind(H(i, j))+0.0121;
+
+% % % % % % % 理想化的Saastamoinen模型/对流层改正
+%             P0 = 1013.25; e0 = 11.69; T0 = 288.15;
+% %           求大气压强
+%             p = P0 * (1.0 - 0.0068 / T0 / apH)^5;
+% %           求开尔文温度
+%             T = 288.15 - 0.0068 * apH;
+% %           水汽压计算
+%             if apH  > 11000
+%                 e = 0;
+%             end
+%             if apH <= 11000
+%                  e = e0 * (1.0 - 0.0068 * apH / T0)^4;
+%             end
+% %           静力学延迟
+%             dtrop = 0.002277 * p / (1.0 - 0.00266 * cosd(2 * H(i, j)) - 0.00028 * apH * 1e-3) ;
+% % % % % % % % % % % % % % % % % % % % % % % % %
+
+
+            if P_t2(i, j) ~= 0
+%           双频组合消除电离层延迟
+            P_t = (f1^2) / (f1^2 - f2^2) * P_t1(i, j) - (f2^2) / (f1^2 - f2^2) * P_t2(i, j);
+            end
+            if P_t2(i, j) == 0
+            P_t = P_t1(i, j);
+            end
+            l(j, 1) = P_t - R + C_V * Dt(i,j) - dtrop;
+            if H(i, j)<5;
+               H(i, j) = H(i, j);
+            end;
+
+            P(j, j) = sind(H(i, j))^2;
+        end
+
+        apX_R = apX;%迭代坐标传参
+        apY_R = apY;
+        apZ_R = apZ;
+%       平差计算得到改正数并进行迭代计算
+        while abs(apX_R - X_R(1,1)) > 0.0001 && abs(apY_R - X_R(2,1)) > 0.0001 && abs(apZ_R - X_R(3,1)) > 0.0001
+            apX_R = X(1,1);
+            apY_R = X(2,1);
+            apZ_R = X(3,1);
+            Q=inv(B.'*P*B);
+            D_X=Q*B.'*P*l;%坐标改正数
+            X(1,1) = apX + D_X(1,1);
+            X(2,1) = apY + D_X(2,1);
+            X(3,1) = apZ + D_X(3,1);%改正后测站坐标
+            X_R(1,1) = X(1,1);
+            X_R(2,1) = X(2,1);
+            X_R(3,1) = X(3,1);
+             
+        end
+        GDOP(1,i) = sqrt(trace(Q));
+
+        StationCLK(1,i) = D_X(4, 1);
+        DeltaXYZ(1,1) = X(1,1) - apX;
+        DeltaXYZ(2,1) = X(2,1) - apY;
+        DeltaXYZ(3,1) = X(3,1) - apZ;
+
+        S = [
+                      -sin(apL),           cos(apL),        0;
+             -sin(apB)*cos(apL), -sin(apB)*sin(apL), cos(apB);
+              cos(apB)*cos(apL),  cos(apB)*sin(apL), sin(apB) 
+            ];
+        
+        ENU_R = S * DeltaXYZ;
+        ENU(1,i) = ENU_R(1,1);
+        ENU(2,i) = ENU_R(2,1);
+        ENU(3,i) = ENU_R(3,1);
+      
+    end
+end
+

Observation_Station_Site-Solving_by_Matlab/Station_Site_Solving.m

@@ -0,0 +1,170 @@
+% 打开坐标文件
+file=fopen('LLA_result_for_Matlab.txt');       
+sPRN = [];
+satX = [];
+satY = [];
+satZ = [];
+P_t1 = [];
+P_t2 = [];
+H    = [];
+TGD  = [];
+Dt   = [];
+sat_num = [];
+
+C_V    = 299792458;
+f1     =   1575.42;
+f2     =   1227.60;
+satnum =         0;
+eponum =         0;
+
+% 设置测站粗略坐标和精确的经纬高数据
+apX = -4052052.8310;
+apY =  4212835.9820;
+apZ = -2545104.4760;
+
+apB = -0.4131213;%纬度
+apL =  2.3367431;%经度
+apH =    603.235;
+
+% 逐行读取文件
+while ~feof(file)
+        line = fgets(file);
+        flag_epo = extractBetween(line,1 ,1 );
+%   判断历元标识
+    if strcmp(flag_epo, ">") == 1
+         eponum = str2double(extractBetween(line,2 ,5 ));
+       while eponum <= 2880 
+               flag_sat = extractBetween(line,1 ,1);
+%          判断卫星标识，读取卫星PRN号和XYZ坐标
+           if strcmp(flag_sat, "G") == 1
+                satnum = satnum + 1;
+                P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16               ,25 + 16 + 16 + 16 ));
+                if P_t1(eponum, satnum) ~= 0
+                    sPRN(eponum, satnum) = str2double(extractBetween(line,2 ,3 ));
+                    satX(eponum, satnum) = str2double(extractBetween(line,11                                   ,25 ));
+                    satY(eponum, satnum) = str2double(extractBetween(line,11 + 16                              ,25 + 16 ));
+                    satZ(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16                         ,25 + 16 + 16 ));
+                    P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16                    ,25 + 16 + 16 + 16 ));
+                    P_t2(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16               ,25 + 16 + 16 + 16 + 16));
+                    H(eponum, satnum)    = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16          ,25 + 16 + 16 + 16 + 16 + 16));
+                    TGD(eponum, satnum)  = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16     ,25 + 16 + 16 + 16 + 16 + 16 + 16));
+                    Dt(eponum, satnum)   = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16 + 16,25 + 16 + 16 + 16 + 16 + 16 + 16 + 16));
+                    line = fgets(file);
+                end
+                if P_t1(eponum, satnum) == 0
+                    satnum = satnum - 1;
+                    line = fgets(file);
+                end
+           end
+%          开头字符不是G，但是是>,判断为一个历元的起点，进行数据读取
+           if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 1
+                satnum = 0;
+                line = fgets(file);
+                continue;
+           end
+%          开头字符既不是G，也不是>,判断为历元结束标志，跳出循环
+           if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 0 && ~feof(file)
+                sat_num(eponum) = satnum;
+                satnum = 0;
+                break;
+           end
+       end
+    end
+    if strcmp(flag_epo, ">") == 0
+        continue;
+    end
+end
+fclose(file);
+
+% 平差计算测站位置和接收机钟差解算
+for i = 1:1:eponum
+    if sat_num(i) >= 4 
+        X = zeros(4, 1);
+        P = zeros(sat_num(i), sat_num(i));
+        B = zeros(sat_num(i), 4);
+        l = zeros(sat_num(i), 1);
+        X_R = zeros(3, 1);
+        
+        DeltaXYZ = zeros(3, 1);
+%       每个历元构成系数矩阵
+        for j = 1:1:sat_num(i)
+            R = sqrt((satX(i, j)-apX)^2+(satY(i, j)-apY)^2+(satZ(i, j)-apZ)^2);
+            B(j, 1) = -((satX(i, j) - apX) / R);
+            B(j, 2) = -((satY(i, j) - apY) / R);
+            B(j, 3) = -((satZ(i, j) - apZ) / R);
+            B(j, 4) = 1;
+
+            dtrop = 2.47/sind(H(i, j))+0.0121;
+
+% % % % % % % 理想化的Saastamoinen模型/对流层改正
+%             P0 = 1013.25; e0 = 11.69; T0 = 288.15;
+% %           求大气压强
+%             p = P0 * (1.0 - 0.0068 / T0 / apH)^5;
+% %           求开尔文温度
+%             T = 288.15 - 0.0068 * apH;
+% %           水汽压计算
+%             if apH  > 11000
+%                 e = 0;
+%             end
+%             if apH <= 11000
+%                  e = e0 * (1.0 - 0.0068 * apH / T0)^4;
+%             end
+% %           静力学延迟
+%             dtrop = 0.002277 * p / (1.0 - 0.00266 * cosd(2 * H(i, j)) - 0.00028 * apH * 1e-3) ;
+% % % % % % % % % % % % % % % % % % % % % % % % %
+
+
+            if P_t2(i, j) ~= 0
+%           双频组合消除电离层延迟
+            P_t = (f1^2) / (f1^2 - f2^2) * P_t1(i, j) - (f2^2) / (f1^2 - f2^2) * P_t2(i, j);
+            end
+            if P_t2(i, j) == 0
+            P_t = P_t1(i, j);
+            end
+            l(j, 1) = P_t - R + C_V * Dt(i,j) - dtrop;
+            if H(i, j)<5;
+               H(i, j) = H(i, j);
+            end;
+
+            P(j, j) = sind(H(i, j))^2;
+        end
+
+        apX_R = apX;%迭代坐标传参
+        apY_R = apY;
+        apZ_R = apZ;
+%       平差计算得到改正数并进行迭代计算
+        while abs(apX_R - X_R(1,1)) > 0.0001 && abs(apY_R - X_R(2,1)) > 0.0001 && abs(apZ_R - X_R(3,1)) > 0.0001
+            apX_R = X(1,1);
+            apY_R = X(2,1);
+            apZ_R = X(3,1);
+            Q=inv(B.'*P*B);
+            D_X=Q*B.'*P*l;%坐标改正数
+            X(1,1) = apX + D_X(1,1);
+            X(2,1) = apY + D_X(2,1);
+            X(3,1) = apZ + D_X(3,1);%改正后测站坐标
+            X_R(1,1) = X(1,1);
+            X_R(2,1) = X(2,1);
+            X_R(3,1) = X(3,1);
+             
+        end
+        GDOP(1,i) = sqrt(trace(Q));
+
+        StationCLK(1,i) = D_X(4, 1);
+        DeltaXYZ(1,1) = X(1,1) - apX;
+        DeltaXYZ(2,1) = X(2,1) - apY;
+        DeltaXYZ(3,1) = X(3,1) - apZ;
+
+        S = [
+                      -sin(apL),           cos(apL),        0;
+             -sin(apB)*cos(apL), -sin(apB)*sin(apL), cos(apB);
+              cos(apB)*cos(apL),  cos(apB)*sin(apL), sin(apB) 
+            ];
+        
+        ENU_R = S * DeltaXYZ;
+        ENU(1,i) = ENU_R(1,1);
+        ENU(2,i) = ENU_R(2,1);
+        ENU(3,i) = ENU_R(3,1);
+      
+    end
+end
+

cal_mod.h

@@ -0,0 +1,271 @@
+#ifndef _CAL_MOD_H_
+#define _CAL_MOD_H_
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+
+#include"read.h"
+#include"2gpst.h"
+#include"select_epo.h"
+
+FILE * result_file = NULL;
+
+#define a  6378137.0//长半轴
+#define f  (1 / 298.257223563)//扁率
+#define e2  (f*(2-f))//第一偏心率平方
+#define C_V 299792458//光速（m）
+#define GM 398600500000000//地心引力常数
+#define math_e 2.718281828459 //e值
+#define PI 3.141592653589793
+#define Earth_e 7.2921151467e-5 //地球自转角速度
+#define C1 0
+
+    //卫星位置结构体
+    typedef struct
+    {
+        double X[36];
+        double Y[36];
+        double Z[36];
+        double deltat[36];//改正前的信号传播时间
+        double delta_t[36];//改正后的信号传播时间
+        double delta_clk[36];
+    }pos_t, *ppos_t;
+
+    //测站位置结构体
+    typedef struct
+    {
+        double X;
+        double Y;
+        double Z;
+        double B;
+        double L;
+        double H;
+        double delta_TR;
+    }station, *pstation;
+
+    typedef struct
+    {
+        double B;
+        int B_d;
+        int B_m;
+        int B_s;
+        double L;
+        int L_d;
+        int L_m;
+        int L_s;
+        double H;
+    }blh, *pblh;
+   
+    typedef struct enu
+    {
+	    double E;
+	    double N;
+	    double U;
+    }enu, *penu;
+
+    typedef struct rah
+    {
+        double R;
+        double A;
+        double H;
+    }rah, *prah;
+
+    //GPS卫星位置计算函数
+    int gps_pos(int eponum, int sPRN, int best_epoch, double GPSsec, 
+                pnav_body nav_b, pobs_head obs_h, ppos_t pos_t)
+    {    
+        double T_S = GPSsec - pos_t[eponum].deltat[sPRN];//计算信号发射时刻
+        double n_0 = sqrt(GM) / pow(nav_b[best_epoch].sqrtA,3);
+        double n = n_0 + nav_b[best_epoch].deltan; //计算卫星运行平均角速度
+        double tk = T_S - nav_b[best_epoch].TOE;//计算归化时间
+        if (tk > 32400)tk -= 604800;
+	    else if (tk < -32400)tk += 604800;//规划时间改正
+        double Ms = nav_b[best_epoch].M0 + n * tk;//计算卫星发射时卫星的平近角M
+	    double Es = Ms, E0;//迭代计算偏近点角
+        do
+	    {
+		E0 = Es;
+		Es = Ms + nav_b[best_epoch].e * sin(Es);
+	    } while (fabs(Es - E0) > 1.0e-5);//迭代计算偏近点角
+        //计算真近点角
+        double cosfs = (cos(Es) - nav_b[best_epoch].e) / (1 - nav_b[best_epoch].e * cos(Es));
+        double sinfs = (sqrt(1 - pow(nav_b[best_epoch].e, 2)) * sin(Es)) / (1 - nav_b[best_epoch].e * cos(Es));
+        double fs = atan2(sinfs, cosfs);
+        double u0 = fs + nav_b[best_epoch].omega;//计算升交角距
+	    //计算摄动改正项
+        double epsilon_u = nav_b[best_epoch].Cus * sin(2 * u0) + nav_b[best_epoch].Cuc * cos(2 * u0);
+        double epsilon_r = nav_b[best_epoch].Crs * sin(2 * u0) + nav_b[best_epoch].Crc * cos(2 * u0);
+        double epsilon_i = nav_b[best_epoch].Cis * sin(2 * u0) + nav_b[best_epoch].Cic * cos(2 * u0);
+        double u = u0 + epsilon_u;//改正后的升交角距
+        double r = pow(nav_b[best_epoch].sqrtA, 2) * (1 - nav_b[best_epoch].e * cos(Es)) + epsilon_r;//改正后的轨道向径
+        double i = nav_b[best_epoch].i0 + epsilon_i + nav_b[best_epoch].IDOT * tk;//改正后的轨道倾角
+        double x = r * cos(u);
+        double y = r * sin(u);//卫星在轨道平面坐标系中的坐标
+        double l = nav_b[best_epoch].OMEGA + (nav_b[best_epoch].deltaomega - Earth_e) * tk - Earth_e * nav_b[best_epoch].TOE;
+        //计算升交点经度
+        double X = x * cos(l) - y * cos(i) * sin(l);
+        double Y = x * sin(l) + y * cos(i) * cos(l);
+        double Z = y * sin(i);//卫星在地固系下的直角坐标
+        pos_t[eponum].X[sPRN] = cos(Earth_e * pos_t[eponum].delta_t[sPRN]) * X + sin(Earth_e * pos_t[eponum].delta_t[sPRN]) * Y;
+        pos_t[eponum].Y[sPRN] = -sin(Earth_e * pos_t[eponum].delta_t[sPRN]) * X + cos(Earth_e * pos_t[eponum].delta_t[sPRN]) * Y;
+        pos_t[eponum].Z[sPRN] = Z;
+        double R = sqrt(pow(pos_t[eponum].X[sPRN] - obs_h->apX, 2) + pow(pos_t[eponum].Y[sPRN] - obs_h->apY, 2) + pow(pos_t[eponum].Z[sPRN] - obs_h->apY, 2));
+        pos_t[eponum].deltat[sPRN] = pos_t[eponum].delta_t[sPRN];
+        pos_t[eponum].delta_t[sPRN] = R / C_V;
+        double rela = 2*sqrt(GM) * nav_b[best_epoch].e * nav_b[best_epoch].sqrtA * sin(Es)/pow(C_V,2);
+        pos_t[eponum].delta_clk[sPRN] = nav_b[best_epoch].sa0 + nav_b[best_epoch].sa1 * tk + nav_b[best_epoch].sa2 * pow(tk,2) - rela;
+    }
+
+    //GPS卫星位置计算
+    int sat_gps_pos_clac(pnav_body nav_b, pobs_epoch obs_e, pobs_body obs_b,
+                         pobs_head obs_h, pstation station, ppos_t pos_t, 
+                         pblh blh, penu enu, prah rah, int o_epochnum, int gps_satnum)
+    {
+        for (int i = 0; i < o_epochnum; i++)//第i个历元
+        {
+            result_file = fopen(".\\Pos_out\\LLA_result_for_read.txt", "a+");
+            fprintf(result_file, "\r>%04d %04d %02d %02d %02d %02d %07.04f",i + 1 ,obs_e[i].y ,obs_e[i].m ,obs_e[i].d ,obs_e[i].h ,obs_e[i].min ,obs_e[i].sec);
+            fclose(result_file);
+            result_file = fopen(".\\Draw the Track[by PyGMT]\\LLA_result_for_PyGMT.txt", "a+");
+            fprintf(result_file, "\r>%04d %04d %02d %02d %02d %02d %07.04f",i + 1 ,obs_e[i].y ,obs_e[i].m ,obs_e[i].d ,obs_e[i].h ,obs_e[i].min ,obs_e[i].sec);
+            fclose(result_file);
+            result_file = fopen(".\\Observation_Station_Site-Solving_by_Matlab\\LLA_result_for_Matlab.txt", "a+");
+            fprintf(result_file, "\n>%04d %04d %02d %02d %02d %02d %07.04f",i + 1 ,obs_e[i].y ,obs_e[i].m ,obs_e[i].d ,obs_e[i].h ,obs_e[i].min ,obs_e[i].sec);
+            fclose(result_file);
+
+            for (int j = 0; j < obs_e[i].gps_num; j++)//第j颗GPS卫星
+            {
+                int y = obs_e[i].y; int m = obs_e[i].m; int d = obs_e[i].d; 
+                double h = obs_e[i].h; int min = obs_e[i].min; double sec = obs_e[i].sec;
+                //转换为GPS周内秒
+                double GPSsec = TimetoGPSsec(y, m, d, h, min, sec);
+                int sPRN = obs_e[i].sPRN_GPS[j];
+                //遍历N文件GPS卫星数据块，寻找最佳历元
+                int best_epoch = select_epoch(GPSsec, sPRN, nav_b, gps_satnum);
+                //观测时刻 - 参考时刻
+                double detat_toc = GPSsec - nav_b[best_epoch].TOE;
+                station[i].delta_TR == 0.0;
+                
+                //计算近似的信号传播时间,接收机钟差已初始化为0(伪距/光速-接收机钟差+卫星钟差)
+                pos_t[i].delta_t[sPRN] = (obs_b[i].obs_gps[j][C1] / C_V) - station[i].delta_TR + nav_b[best_epoch].sa0 + nav_b[best_epoch].sa1 * detat_toc + nav_b[best_epoch].sa2 * pow(detat_toc, 2);
+                
+                //判断收敛
+                pos_t[i].deltat[sPRN] = 0.0;
+                while (fabs(pos_t[i].delta_t[sPRN] - pos_t[i].deltat[sPRN]) > 0.00001){
+                    
+                    gps_pos( i, sPRN, best_epoch, GPSsec, nav_b, obs_h, pos_t);
+                }
+
+                //地固坐标系转经纬大地高坐标系******************************************************//
+                double X = pos_t[i].X[sPRN];
+                double Y = pos_t[i].Y[sPRN];
+                double Z = pos_t[i].Z[sPRN];
+                double B = 0.0, N = 0.0, H = 0.0, R0, R1, deltaH, deltaB;
+                R0 = sqrt(pow(X, 2) + pow(Y, 2));
+                R1 = sqrt(pow(X, 2) + pow(Y, 2) + pow(Z, 2));
+                //经度直接求解
+                blh->L = atan2(Y, X) * 57.295779513;
+                //迭代求大地维度和大地高
+                N = a;
+                H = R1 - a;
+                B = atan2(Z * (N + H), R0 * (N * (1 - e2) + H));
+                do
+                {
+                    deltaH = N;//判断收敛所用
+                    deltaB = B;
+		            N = a / sqrt(1 - e2 * pow(sin(B), 2));
+		            H = R0 / cos(B) - N;
+		            B = atan2(Z * (N + H), R0 * (N * (1 - e2) + H));
+	            } while (fabs(deltaH - H) > 1.0e-3 && fabs(deltaB - B) > 1.0e-9);
+	            blh->B = B * 57.295779513;
+	            blh->H = H ;
+                //*******************************************************************************//
+
+                //地固坐标系转经纬大地高坐标系******************************************************//
+                X = obs_h->apX;
+                Y = obs_h->apY;
+                Z = obs_h->apZ;
+                B = 0.0, N = 0.0, H = 0.0;
+                R0 = sqrt(pow(X, 2) + pow(Y, 2));
+                R1 = sqrt(pow(X, 2) + pow(Y, 2) + pow(Z, 2));
+                //经度直接求解
+                station->L = atan2(Y, X);
+                //迭代求大地维度和大地高
+                N = a;
+                H = R1 - a;
+                B = atan2(Z * (N + H), R0 * (N * (1 - e2) + H));
+                do
+                {
+                    deltaH = N;//判断收敛所用
+                    deltaB = B;
+		            N = a / sqrt(1 - e2 * pow(sin(B), 2));
+		            H = R0 / cos(B) - N;
+		            B = atan2(Z * (N + H), R0 * (N * (1 - e2) + H));
+	            } while (fabs(deltaH - H) > 1.0e-3 && fabs(deltaB - B) > 1.0e-9);
+	            station->B = B;
+	            station->H = H;
+                //*******************************************************************************//
+
+                //大地高坐标系转东北天坐标系*******************************************************//
+                double sinL = sin(station->L);
+                double cosL = cos(station->L);
+                double sinB = sin(station->B);
+                double cosB = cos(station->B);
+                enu->E = -sinL*(pos_t[i].X[sPRN] - obs_h->apX) + cosL*(pos_t[i].Y[sPRN] - obs_h->apY);
+                enu->N = -sinB*cosL*(pos_t[i].X[sPRN] - obs_h->apX) - sinB*sinL*(pos_t[i].Y[sPRN] - obs_h->apY) + cosB*(pos_t[i].Z[sPRN] - obs_h->apZ);
+                enu->U = cosB*cosL*(pos_t[i].X[sPRN] - obs_h->apX) + cosB*sinL*(pos_t[i].Y[sPRN] - obs_h->apY) + sinB*(pos_t[i].Z[sPRN] - obs_h->apZ);
+                //*******************************************************************************//
+
+                //卫星方位角a，高度角h，向径r计算**************************************************//
+                rah->H = atan2(enu->U, sqrt(enu->E * enu->E + enu->N * enu->N)) * 57.295779513;
+                rah->A = atan2(enu->E, enu->U) * 57.295779513;
+                    if (rah->A < 0)
+                        rah->A += 2 * PI;
+                    if (rah->A > 2 * PI)
+                        rah->A -= 2 * PI;
+                rah->R = sqrt(enu->E * enu->E + enu->N * enu->N + enu->U * enu->U);
+
+                //*******************************************************************************//
+                if (blh->H < 0 || nav_b[best_epoch].sHEA != 0)
+                { 
+                    break;
+                }
+                else
+                {  int B_sig = blh->B / fabs(blh->B);//符号判断
+                    blh->B_d = (int)(blh->B);//转换为度分秒
+                    blh->B_m = (int)((blh->B - blh->B_d) * 60);//分秒去除符号格式
+                    blh->B_s = (int)(((blh->B - blh->B_d) * 60 - blh->B_m) * 60);
+                    blh->B_m = blh->B_m * B_sig;
+                    blh->B_s = blh->B_s * B_sig;
+
+                    int L_sig = blh->L / fabs(blh->L);
+                    blh->L_d = (int)(blh->L);//转换为度分秒
+                    blh->L_m = (int)((blh->L - blh->L_d) * 60);
+                    blh->L_s = (int)(((blh->L - blh->L_d) * 60 - blh->L_m) * 60);
+                    blh->L_m = blh->L_m * L_sig;
+                    blh->L_s = blh->L_s * L_sig;
+
+                    result_file = fopen(".\\Pos_out\\LLA_result_for_read.txt", "a+");
+                    fprintf(result_file, "\nG%02d[sPRN] %4d° %02d′ %02d″ %4d° %02d′ %02d″ %15.05f",sPRN, blh->B_d, blh->B_m, blh->B_s, blh->L_d, blh->L_m, blh->L_s, blh->H);
+                    fclose(result_file);
+                    result_file = fopen(".\\Draw the Track[by PyGMT]\\LLA_result_for_PyGMT.txt", "a+");
+                    fprintf(result_file, "\nG%02d[sPRN] %10.05f %10.05f %15.05f",sPRN, blh->B, blh->L, blh->H);
+                    fclose(result_file);
+                    result_file = fopen(".\\Observation_Station_Site-Solving_by_Matlab\\LLA_result_for_Matlab.txt", "a+");
+                    fprintf(result_file, "\nG%02d[sPRN] %15.05f %15.05f %15.05f %15.05f %15.05f %15.05f %15.12f %16.13f"
+                                            ,sPRN, pos_t[i].X[sPRN], pos_t[i].Y[sPRN], pos_t[i].Z[sPRN], obs_b[i].obs_gps[j][C1], obs_b[i].obs_gps[j][10], rah->H, nav_b[best_epoch].TGD, pos_t[i].delta_clk[sPRN]);
+                    fclose(result_file);
+                }
+            }       
+        
+            result_file = fopen(".\\Observation_Station_Site-Solving_by_Matlab\\LLA_result_for_Matlab.txt", "a+");
+            fprintf(result_file, "\n");
+            fclose(result_file);
+        }
+        
+        result_file = fopen(".\\Observation_Station_Site-Solving_by_Matlab\\LLA_result_for_Matlab.txt", "a+");
+        fprintf(result_file, "\nEND");
+        fclose(result_file);
+    }
+
+#endif

main.c

@@ -0,0 +1,64 @@
+//仅支持{ Ver 3.04 }版本RENIX观测值文件与广播星历文件读取
+
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+
+#include"read.h"
+#include"cal_mod.h"
+
+int main(void)
+{
+	//数据读取
+	FILE * fp_nav = NULL;//导航星历文件指针
+	FILE * fp_obs = NULL;//观测值文件指针
+
+	pstation station = NULL;
+	ppos_t pos_t = NULL;
+	pblh blh = NULL;
+	penu enu = NULL;
+	prah rah = NULL;
+
+	pnav_head nav_h = NULL;
+	pnav_body nav_b = NULL;
+	pobs_head obs_h = NULL;
+	pobs_epoch obs_e = NULL;
+	pobs_body obs_b = NULL;//创建结构体指针，并将其初始化为空
+
+	//O文件读取
+	fp_obs = fopen("alic0600.22o", "r");//以只读的方式打开O文件
+	int o_epochnum = get_epochnum(fp_obs);//O文件历元数
+	rewind(fp_obs);
+	obs_h = (pobs_head)malloc(sizeof(obs_head));//给O文件头开辟空间
+	obs_e = (pobs_epoch)malloc(sizeof(obs_epoch) * o_epochnum);
+	obs_b = (pobs_body)malloc(sizeof(obs_body) * o_epochnum);
+		if (obs_h && obs_e && obs_b)
+		{
+			read_o_h(fp_obs, obs_h);//读取O文件头
+			read_o_b(fp_obs, obs_e, obs_b, obs_h->obstypenum_gps, obs_h->obstypenum_bds);//读取O数据块
+		}
+		fclose(fp_obs);//关闭O文件
+
+	//N文件读取
+	fp_nav = fopen("brdm0600.22p", "r");//以只读的方式打开N文件
+	int gps_satnum = getgpssatnum(fp_nav);//获取导航文件GPS卫星数据块数
+	rewind(fp_nav);//将文件指针返回值起始位置
+	nav_h = (pnav_head)malloc(sizeof(nav_head));//给N文件头开辟空间
+	nav_b = (pnav_body)malloc(sizeof(nav_body) * gps_satnum * 20);
+		if (nav_h && nav_b)
+		{
+			read_n_h(fp_nav, nav_h);//读取N头文件
+			read_n_b(fp_nav, nav_b);//读取N数据块
+		}
+		fclose(fp_nav);//关闭N文件
+
+	blh = (pblh)malloc(sizeof(pblh) * o_epochnum *256);
+	enu = (penu)malloc(sizeof(penu) * o_epochnum *256);
+	rah = (prah)malloc(sizeof(prah) * o_epochnum *256);
+	pos_t = (ppos_t)malloc(sizeof(ppos_t) * o_epochnum * 256);
+	station = (pstation)malloc(sizeof(pstation) * o_epochnum * 256);
+	sat_gps_pos_clac(nav_b, obs_e, obs_b, obs_h, station, pos_t, blh, enu, rah, o_epochnum, gps_satnum);
+
+	
+	return 0;
+}

read.h

@@ -0,0 +1,601 @@
+#ifndef _READ_H_
+#define _READ_H_
+
+#include<stdio.h>
+#include<string.h>
+
+#define MAXRINEX 1000//一行最大字符数
+
+void setstr(char* des, const char* src, int n)
+{
+	char* p = des;
+	const char* q = src;
+	while (*q && q < src+n)
+	{
+		*p++ = *q++;
+	}
+	*p-- = '\0';
+//去掉尾部空格
+	while (p >= des && *p == ' ')
+	{
+		*p-- = '\0';
+	}
+}
+
+//获取文件数据块行数，从END OF HEADER后开始起算
+extern int getgpssatnum(FILE* fp_nav)
+{
+	int gps_satnum = 0;
+	int flag = 0;
+	char buff[MAXRINEX];//用来存放读取到的字符串
+	char satvar;
+	char* lable = buff + 60;
+	//fgets函数，读取一行，当读取结束后返回NULL指针，格式如下：
+	//char * fgets ( char * str, int num, FILE * stream );
+	while (fgets(buff, MAXRINEX, fp_nav))
+	{
+		if (flag == 1)
+		{
+			while (fgets(buff, MAXRINEX, fp_nav))
+			{
+				strncpy(&satvar, buff + 0, 1);
+				if (satvar != 'G')
+				{
+					continue;
+				}
+				else if (satvar == 'G')
+				{
+					gps_satnum++;
+					break;	
+				}	
+			}
+											
+		}
+
+		if (strstr(lable, "END OF HEADER"))
+		{
+			flag = 1;
+		}
+	}
+	return gps_satnum;
+}
+
+//将字符串转换为浮点数,i起始位置，n输入多少个字符
+static double strtonum(const char* buff, int i, int n)
+{
+	double value = 0.0;
+	char str[256] = { 0 };
+	char* p = str;
+	/************************************
+	* 当出现以下三种情况报错，返回0.0
+	* 1.起始位置<0
+	* 2.读取字符串个数<i
+	* 3.str里面存放的字节数<n
+	*************************************/
+	if (i < 0 || (int)strlen(buff) < i || (int)sizeof(str) - 1 < n)
+	{
+		return 0.0;
+	}
+	for (buff += i; *buff && --n >= 0; buff++)
+	{
+		//三目操作符：D和d为文件中科学计数法部分，将其转换成二进制能读懂的e
+		*p++ = ((*buff == 'D' || *buff == 'd') ? 'e' : *buff);
+	}
+	*p = '\0';
+	//三目操作符，将str中存放的数以格式化读取到value中。
+	return sscanf(str, "%lf", &value) == 1 ? value : 0.0;
+}
+
+//导航头文件结构体
+typedef struct nav_head
+{
+	double ver;//rinex 版本号
+	char type[20];//读取的数据类型
+	double ION_GPSA[4];//8个电离层参数
+	double ION_GPSB[4];
+	double ION_BDSA[4];//8个电离层参数
+	double ION_BDSB[4];
+	double ION_QZSA[4];//8个电离层参数
+	double ION_QZSB[4];
+	double ION_IRNA[4];//8个电离层参数
+	double ION_IRNB[4];
+	double GPUT[4];
+	double GLUT[4];
+	double GAUT[4];
+	double BDUT[4];
+	double QZUT[4];
+	double IRUT[4];
+	double SBUT[4];
+	int leap;
+}nav_head, *pnav_head;
+
+//导航数据结构体
+typedef struct nav_body
+{
+	//数据块第一行内容：
+	int sPRN;//卫星PRN号
+	//历元：TOC中卫星钟的参考时刻
+	int TOC_Y;//年
+	int TOC_M;//月
+	int TOC_D;//日
+	int TOC_H;//时
+	int TOC_Min;//分
+	int TOC_Sec;//秒
+	double sa0;//卫星钟差
+	double sa1;//卫星钟偏
+	double sa2;//卫星钟漂
+ 
+	//数据块第二行内容：
+	double IODE;//数据、星历发布时间(数据期龄)
+	double Crs;//轨道半径的正弦调和改正项的振幅（单位：m）
+	double deltan;//卫星平均运动速率与计算值之差(rad/s)
+	double M0;//参考时间的平近点角(rad)
+ 
+	//数据块第三行内容：
+	double Cuc;//维度幅角的余弦调和改正项的振幅(rad)
+	double e;//轨道偏心率
+	double Cus;//轨道幅角的正弦调和改正项的振幅(rad)
+	double sqrtA;//长半轴平方根
+ 
+	//数据块第四行内容：
+	double TOE;//星历的参考时刻(GPS周内秒)
+	double Cic;//轨道倾角的余弦调和改正项的振幅(rad)
+	double OMEGA;//参考时刻的升交点赤经
+	double Cis;//维度倾角的正弦调和改正项的振幅(rad)
+ 
+    //数据块第五行内容：
+	double i0;//参考时间的轨道倾角(rad)
+	double Crc;//轨道平径的余弦调和改正项的振幅(m)
+	double omega;//近地点角距
+	double deltaomega;//升交点赤经变化率(rad)
+ 
+	//数据块第六行内容：
+	double IDOT;//近地点角距(rad/s)
+	double L2code;//L2上的码
+	double GPSweek;//GPS周,于TOE一同表示
+	double L2Pflag;//L2,p码数据标记
+ 
+	//数据块第七行内容
+	double sACC;//卫星精度
+	double sHEA;//卫星健康状态
+	double TGD;//sec
+	double IODC;//钟的数据龄期
+ 
+	//数据块第八行内容
+	double TTN;//电文发送时间
+	double fit;//拟合区间
+	double spare1;//空
+	double spare2;//空
+ 
+}nav_body, *pnav_body;
+
+//观测值头文件
+typedef struct obs_head
+{
+	double ver;//RINEX文件版本号
+	char type[30];//文件类型
+	double apX;//测站近似位置XYZ（WGS84）
+	double apY;
+	double apZ;
+	double ANTH;//天线高
+	double ANTdeltaE;//天线中心对于测站标志在东方向上的偏移量
+	double ANTdeltaN;//天线中心对于测站标志在北方向上的偏移量
+	int WAFL1;//缺省的L1载波的波长因子
+	int WAFL2;//缺省的L2载波的波长因子
+	int WAFflag;
+	int obstypenum_gps;//GPS观测值类型数量
+	int obstypenum_bds;
+	char obstype_gps[36];//rinex3.04GPS观测值类型列表
+	char obstype_bds[36];
+	double interval;//观测值的历元间隔
+	int f_y;//第一个观测记录的时刻
+	int f_m;
+	int f_d;
+	int f_h;
+	int f_min;
+	double f_sec;
+	char tsys[5];//时间系统
+}obs_head,*pobs_head;
+
+
+//观测值历元数据结构体
+typedef struct obs_epoch
+{
+	//观测历元时刻
+	int y;
+	int m;
+	int d;
+	int h;
+	int min;
+	double sec;
+	int p_flag;//历元标志
+	int sat_num;//当前历元所观测到的卫星数量
+	int gps_num;//当前历元所观测到的GPS卫星数量
+	int bds_num;//当前历元所观测到的BDS卫星数量
+	int sPRN[84];//当前历元所能观测到的卫星的PRN列表
+	int sPRN_GPS[25];
+	int sPRN_BDS[36];
+}obs_epoch,*pobs_epoch;
+
+
+//观测值数据结构体
+typedef struct obs_body
+{
+	double obs_gps[100][24];//GPS观测值
+	double obs_bds[100][24];//BDS观测值
+}obs_body,*pobs_body;
+
+
+//读取O文件******************************************************************************************************// 
+//获取O文件历元数
+extern int get_epochnum(FILE* fp_obs)
+{
+	int n = 0;//记录历元数
+	int satnum = 0;//记录每个历元的卫星数
+	char flag;//存放标志符号'>'
+	char buff[MAXRINEX];//存放读取的字符串
+	while (fgets(buff, MAXRINEX, fp_obs))
+	{
+		satnum = (int)strtonum(buff, 33, 2);
+		strncpy(&flag, buff + 0, 1);
+		if (flag == '>')
+		{
+			n++;
+		}
+	}
+	return n;
+}
+
+//读取O文件数据头
+extern void read_o_h(FILE* fp_obs, pobs_head obs_h)
+{
+	char buff[MAXRINEX] = { 0 };
+	char flag = { 0 };
+	char* lable = buff + 60;
+	int i = 0;
+	int j = 0;
+	while (fgets(buff, MAXRINEX, fp_obs))
+	{
+		if (strstr(lable, "RINEX VERSION / TYPE"))
+		{
+			obs_h->ver = strtonum(buff, 0, 9);
+			strncpy(obs_h->type, buff + 20, 30);
+			continue;
+		}
+		else if (strstr(lable, "APPROX POSITION XYZ"))
+		{
+			obs_h->apX = strtonum(buff, 0, 14);
+			obs_h->apY = strtonum(buff, 0 + 14, 14);
+			obs_h->apZ = strtonum(buff, 0 + 14 + 14, 14);
+			continue;
+		}
+		else if (strstr(lable, "ANTENNA: DELTA H/E/N"))
+		{
+			obs_h->ANTH = strtonum(buff, 0, 14);
+			obs_h->ANTdeltaE = strtonum(buff, 14, 14);
+			obs_h->ANTdeltaN = strtonum(buff, 14 + 14, 14);
+			continue;
+		}
+		else if (strstr(lable, "WAVELENGTH FACT L1/2"))
+		{
+			obs_h->WAFL1 = (int)strtonum(buff, 0, 6);
+			obs_h->WAFL2 = (int)strtonum(buff, 6, 6);
+			obs_h->WAFflag = (int)strtonum(buff, 6 + 6, 6);
+			continue;
+		}
+		else if (strstr(lable, "SYS / # / OBS TYPES"))
+		{
+			char flag = {0};
+			strncpy(&flag, buff + 0, 1);
+			if (flag == 'G')
+			{
+				obs_h->obstypenum_gps = (int)strtonum(buff, 4, 2);
+				if (obs_h->obstypenum_gps <= 13)
+				{
+					for (i = 0; i < obs_h->obstypenum_gps; i++)
+					{
+						strncpy(&(obs_h->obstype_gps[i]), buff + 7 + 4 * i, 3);
+					}
+				}
+				else if (obs_h->obstypenum_gps > 13)
+				{
+					for (i = 0; i < 13; i++)
+					{
+						strncpy(&(obs_h->obstype_gps[i]), buff + 7 + 4 * i, 3);
+					}
+					fgets(buff, MAXRINEX, fp_obs);
+					for (i = 0; i < obs_h->obstypenum_gps - 13; i++)
+					{
+						strncpy(&(obs_h->obstype_gps[i + 13]), buff + 7 + 4 * i, 3);
+					}
+				}
+				continue;
+			}
+			if (flag == 'C')
+			{
+				obs_h->obstypenum_bds = (int)strtonum(buff, 4, 2);
+				if (obs_h->obstypenum_bds <= 13)
+				{
+					for (i = 0; i < obs_h->obstypenum_bds; i++)
+					{
+						strncpy(&(obs_h->obstype_bds[i]), buff + 7 + 4 * i, 3);
+					}
+				}
+				else if (obs_h->obstypenum_bds > 13)
+				{
+					for (i = 0; i < 13; i++)
+					{
+						strncpy(&(obs_h->obstype_bds[i]), buff + 7 + 4 * i, 3);
+					}
+					fgets(buff, MAXRINEX, fp_obs);
+					for (i = 0; i < obs_h->obstypenum_bds - 13; i++)
+					{
+						strncpy(&(obs_h->obstype_bds[i + 13]), buff + 7 + 4 * i, 3);
+					}
+				}
+				continue;
+			}
+		}
+		else if (strstr(lable, "INTERVAL"))
+		{
+			obs_h->interval = strtonum(buff, 0, 11);
+			continue;
+		}
+		else if (strstr(lable, "TIME OF FIRST OBS"))
+		{
+			obs_h->f_y = (int)strtonum(buff, 2, 4);
+			obs_h->f_m = (int)strtonum(buff, 2 + 6, 4);
+			obs_h->f_d = (int)strtonum(buff, 2 + 6 + 6, 4);
+			obs_h->f_h = (int)strtonum(buff, 2 + 6 + 6 + 6, 4);
+			obs_h->f_min = (int)strtonum(buff, 2 + 6 + 6 + 6 + 6, 4);
+			obs_h->f_sec = strtonum(buff, 2 + 6 + 6 + 6 + 6 + 6, 9);
+			strncpy(obs_h->tsys, buff + 6 + 6 + 6 + 6 + 6 + 18, 3);
+			continue;
+		}
+		else if (strstr(lable, "END OF HEADER"))
+			break;
+	}
+}
+
+//读取O文件数据块
+extern void read_o_b(FILE* fp_obs, pobs_epoch obs_e, pobs_body obs_b, int type_gps, int type_bds)
+{
+	int n = 0;//历元数
+	char buff[MAXRINEX] = { 0 };
+	char flag = { 0 };//判断符号
+	while (fgets(buff, MAXRINEX, fp_obs))
+	{
+		int i = 0;//第i颗卫星
+	    int j_g = 0;//第j颗GPS卫星
+		int k_g = 0;//第j颗GPS卫星的第k个观测值
+		int j_c = 0;//第j颗BDS卫星
+		int k_c = 0;//第j颗BDS卫星的第k个观测值
+		obs_e[n].gps_num = 0;
+		obs_e[n].bds_num = 0;
+
+		//按照格式将历元参考时间依次存入
+		//除秒以外，均转换成整型
+		obs_e[n].y = (int)strtonum(buff, 2, 4);
+		obs_e[n].m = (int)strtonum(buff, 7, 2);
+		obs_e[n].d = (int)strtonum(buff, 10, 2);
+		obs_e[n].h = (int)strtonum(buff, 13, 2);
+		obs_e[n].min = (int)strtonum(buff, 16, 2);
+		obs_e[n].sec = strtonum(buff, 19, 10);
+		obs_e[n].p_flag = (int)strtonum(buff, 31, 1);
+		//输入卫星数
+		obs_e[n].sat_num = strtonum(buff, 33, 3);
+		strncpy(&flag, buff + 0, 1);
+		
+		for (i = 0; i < obs_e[n].sat_num; i++)
+		{									
+			//判断是否为GPS卫星
+			if (flag == 'G')//获取flag标识，判断GPS卫星为真，读取该卫星PRN号
+			{		
+				obs_e[n].sPRN[i] = strtonum(buff, 1, 2);//将卫星编号存储到第n个历元中，第i个PRN位置
+				obs_e[n].sPRN_GPS[j_g] = obs_e[n].sPRN[i];//将GPS卫星编号与总观测卫星号关联
+				for (k_g = 0; k_g < type_gps; k_g++)
+				{
+					obs_b[n].obs_gps[j_g][k_g] = strtonum(buff, 3 + 16 * k_g, 14);
+				}				
+				obs_e[n].gps_num++; j_g++;
+				fgets(buff, MAXRINEX, fp_obs);
+				strncpy(&flag, buff + 0, 1);				
+			}
+						
+			else if (flag != 'G')//获取flag标识，判断GPS卫星为否
+			{	 
+				//判断是否为BDS卫星
+				 if (flag == 'C')//获取flag标识，判断BDS卫星为真，读取该卫星PRN号
+				{
+					obs_e[n].sPRN[i] = strtonum(buff, 1, 2);//将卫星编号存储到第n个历元中，第i个PRN位置
+				 	obs_e[n].sPRN_BDS[j_c] = obs_e[n].sPRN[i];//将BDS卫星编号与总观测卫星号关联
+					for (k_c = 0; k_c < type_bds; k_c++)
+					{
+						obs_b[n].obs_bds[j_c][k_c] = strtonum(buff, 3 + 16 * k_c, 14);
+					}				 
+				 	obs_e[n].bds_num++; j_c++;
+				 	fgets(buff, MAXRINEX, fp_obs);
+				 	strncpy(&flag, buff + 0, 1);
+				}				
+				 
+				else if (flag != 'C')//获取flag标识，判断BDS卫星为真，读取该卫星PRN号
+				{
+				 	fgets(buff, MAXRINEX, fp_obs);
+					strncpy(&flag, buff + 0, 1);
+				}
+			}
+
+		}		
+	n++;
+	}
+	
+}
+
+//读取N文件***************************************************************************************************************//
+//读取N文件数据头
+void read_n_h(FILE* fp_nav, pnav_head nav_h)
+{
+	char buff[MAXRINEX] = { 0 };
+	char* lable = buff;
+	int i = 0;
+	int j = 0;
+	while (fgets(buff, MAXRINEX, fp_nav))
+	{
+		if (strstr(lable, "RINEX VERSION / TYPE"))
+		{
+			nav_h->ver = strtonum(buff, 0, 9);
+			strncpy((nav_h->type), buff + 20, 15);
+			continue;
+		}
+		else if (strstr(lable, "GPSA"))
+		{
+			nav_h->ION_GPSA[0] = strtonum(buff, 6, 12);
+			nav_h->ION_GPSA[1] = strtonum(buff, 6 + 12, 12);
+			nav_h->ION_GPSA[2] = strtonum(buff, 6 + 12 + 12, 12);
+			nav_h->ION_GPSA[3] = strtonum(buff, 6 + 12 + 12 + 12, 12);
+			continue;
+		}
+		else if (strstr(lable, "GPSB"))
+		{
+			nav_h->ION_GPSB[0] = strtonum(buff, 6, 12);
+			nav_h->ION_GPSB[1] = strtonum(buff, 6 + 12, 12);
+			nav_h->ION_GPSB[2] = strtonum(buff, 6 + 12 + 12, 12);
+			nav_h->ION_GPSB[3] = strtonum(buff, 6 + 12 + 12 + 12, 12);
+			continue;
+		}
+		else if (strstr(lable, "BDSA"))
+		{
+			nav_h->ION_BDSA[0] = strtonum(buff, 6, 12);
+			nav_h->ION_BDSA[1] = strtonum(buff, 6 + 12, 12);
+			nav_h->ION_BDSA[2] = strtonum(buff, 6 + 12 + 12, 12);
+			nav_h->ION_BDSA[3] = strtonum(buff, 6 + 12 + 12 + 12, 12);
+			continue;
+		}
+		else if (strstr(lable, "BDSB"))
+		{
+			nav_h->ION_BDSB[0] = strtonum(buff, 6, 12);
+			nav_h->ION_BDSB[1] = strtonum(buff, 6 + 12, 12);
+			nav_h->ION_BDSB[2] = strtonum(buff, 6 + 12 + 12, 12);
+			nav_h->ION_BDSB[3] = strtonum(buff, 6 + 12 + 12 + 12, 12);
+			continue;
+		}
+		else if (strstr(lable, "BDUT"))
+		{
+			nav_h->BDUT[0] = strtonum(buff, 5, 17);
+			nav_h->BDUT[1] = strtonum(buff, 5 + 17, 16);
+			nav_h->BDUT[2] = strtonum(buff, 5 + 17 + 16, 7);
+			nav_h->BDUT[3] = strtonum(buff, 5 + 17 + 16 + 7, 5);
+			continue;
+		}
+		else if (strstr(lable, "GPUT"))
+		{
+			nav_h->GPUT[0] = strtonum(buff, 5, 17);
+			nav_h->GPUT[1] = strtonum(buff, 5 + 17, 16);
+			nav_h->GPUT[2] = strtonum(buff, 5 + 17 + 16, 7);
+			nav_h->GPUT[3] = strtonum(buff, 5 + 17 + 16 + 7, 5);
+			continue;
+		}
+
+		else if (strstr(lable, "LEAP SECONDS"))
+		{
+			nav_h->leap = (int)strtonum(buff, 4, 2);
+		}
+		else if (strstr(lable, "END OF HEADER"))//头文件读取结束
+		break;
+	}
+}
+
+//读取N文件数据块
+void read_n_b(FILE* fp_nav, pnav_body nav_b)
+{
+	int i_g = 0;//第i_g个GPS卫星数据
+	int i_c = 0;//第i_c个BDS卫星数据
+	char buff[84] = { 0 };
+	char flag = { 0 };//判断符号
+	while (fgets(buff, MAXRINEX, fp_nav))
+	{
+		int j_g = 0;//数据块的第j_g行数据
+		int j_c = 0;//数据块的第j_C行数据
+		strncpy(&flag, buff + 0, 1);
+		if(flag == 'G')
+		{
+			for (j_g = 0; j_g < 8; j_g++)//读取第i_g颗 GPS卫星广播星历("总共包含8行数据") 数据块的第j_g行数据
+			{
+				switch (j_g)
+				{
+				case 0:
+					nav_b[i_g].sPRN = (int)strtonum(buff, 1, 2);
+					nav_b[i_g].TOC_Y = (int)strtonum(buff, 4, 4);
+					nav_b[i_g].TOC_M = (int)strtonum(buff, 9, 2);
+					nav_b[i_g].TOC_D = (int)strtonum(buff, 12, 2);
+					nav_b[i_g].TOC_H = (int)strtonum(buff, 15, 2);
+					nav_b[i_g].TOC_Min = (int)strtonum(buff, 18, 2);
+					nav_b[i_g].TOC_Sec = strtonum(buff, 21, 2);
+					nav_b[i_g].sa0 = strtonum(buff, 23, 19);
+					nav_b[i_g].sa1 = strtonum(buff, 23 + 19, 19);
+					nav_b[i_g].sa2 = strtonum(buff, 23 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 1:
+					nav_b[i_g].IODE = strtonum(buff, 4, 19);
+					nav_b[i_g].Crs = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].deltan = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].M0 = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 2:
+					nav_b[i_g].Cuc = strtonum(buff, 4, 19);
+					nav_b[i_g].e = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].Cus = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].sqrtA = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 3:
+					nav_b[i_g].TOE = strtonum(buff, 4, 19);
+					nav_b[i_g].Cic = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].OMEGA = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].Cis = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 4:
+					nav_b[i_g].i0 = strtonum(buff, 4, 19);
+					nav_b[i_g].Crc = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].omega = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].deltaomega = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 5:
+					nav_b[i_g].IDOT = strtonum(buff, 4, 19);
+					nav_b[i_g].L2code = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].GPSweek= strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].L2Pflag = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 6:
+					nav_b[i_g].sACC = strtonum(buff, 4, 19);
+					nav_b[i_g].sHEA = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].TGD = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].IODC = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					fgets(buff, MAXRINEX, fp_nav);
+					break;
+				case 7:
+					nav_b[i_g].TTN = strtonum(buff, 4, 19);
+					nav_b[i_g].fit = strtonum(buff, 4 + 19, 19);
+					nav_b[i_g].spare1 = strtonum(buff, 4 + 19 + 19, 19);
+					nav_b[i_g].spare2 = strtonum(buff, 4 + 19 + 19 + 19, 19);
+					break;
+				}				
+			}
+			i_g++;
+			continue;
+		}
+		else if(flag != 'G')
+		{
+			fgets(buff, MAXRINEX, fp_nav);
+			continue;
+		}		
+	}
+}
+
+
+#endif

select_epo.h

@@ -0,0 +1,31 @@
+#ifndef _SELECT_EPO_H_
+#define _SELECT_EPO_H_
+
+#include<stdio.h>
+#include<math.h>
+
+#include"read.h"
+
+//挑选合适的星历
+extern int select_epoch(double SecofWeek, int sPRN, pnav_body nav_b, int gps_satnum)
+{
+	int best_epoch = 0;
+	double min=10000;//假设最小值是1万秒
+	double Min;
+	int i;
+	for (i = 0; i < gps_satnum; i++)
+	{
+		if (sPRN == nav_b[i].sPRN)
+		{
+			Min = fabs(SecofWeek - nav_b[i].TOE);
+			if (Min <= min)
+			{
+				best_epoch = i;
+				min = Min;
+			}
+		}
+	}
+	return best_epoch;
+}
+
+#endif