package location;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.lang.reflect.Array;
import java.lang.reflect.Constructor;
import java.text.DateFormat;
import java.text.FieldPosition;
import java.text.ParsePosition;
import java.time.DayOfWeek;
import java.time.LocalDateTime;
import java.time.ZoneId;
import java.util.ArrayList;
import java.util.Date;
import java.util.Iterator;
import java.util.ListIterator;
import java.util.Stack;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
public class location {
public static String calculateLocation(String satellite,double sqrt_a,double delta_n,double t,double toe,double M0,double e,double omg,double Cuc,double Cus,double Crc,double Crs,double Cic,double Cis,double i0,double OMG0,double OMG0_DOT) {
//长半径
double a=Math.pow(sqrt_a, 2);
//地球重力常数
double constGM=3.9860047e14;
//平均角速度
double n0=Math.sqrt(constGM/Math.pow(a, 3));
double n=n0+delta_n;
//平近点角
Double tk=t-toe;
double Mk=M0+n*(t-toe);
//偏近点角
double Ek0=0;
double Ek1=Mk;
while(Math.abs(Ek1-Ek0)>10e-12) {
Ek0=Ek1;
Ek1=Mk+e*Math.sin(Ek0);
}
double Ek=Ek0;
//真近点角
double vk=Math.acos((Math.cos(Ek)-e)/(1-e*Math.cos(Ek)));
//升交角距
double uk=vk+omg;
//扰动改正
double delta_uk=Cuc*Math.cos(2*uk)+Cus*Math.sin(2*uk);
double delta_rk=Crc*Math.cos(2*uk)+Crs*Math.sin(2*uk);
double delta_ik=Cic*Math.cos(2*uk)+Cis*Math.sin(2*uk);
//赤经 半径 轨道倾角改正
double u=uk+delta_uk;
double r=a*(1-e*Math.cos(Ek))+delta_rk;
double i=i0+delta_ik;
//卫星在轨道上的位置
double x=r*Math.cos(u);
double y=r*Math.sin(u);
//地球自转角速度
double omge=7.292115e-5;
//历元t升交点的赤纬
double lambda=OMG0+(OMG0_DOT*omge)*tk-omge*toe;
//卫星在地心地固坐标系中的位置
double x_E=r*(Math.cos(u)*Math.cos(lambda)-Math.sin(u)*Math.cos(i)*Math.sin(lambda));
double y_E=r*(Math.cos(u)*Math.sin(lambda)+Math.sin(u)*Math.cos(i)*Math.cos(lambda));
double z_E=r*(Math.sin(u)*Math.sin(i));
return "卫星:"+satellite+",长半径:"+a+",平均角速度:"+n0
+",\n偏近点角:"+Ek+",真近点角"+vk+",升交角据:"+uk
+",\n改正赤经:"+u+",改正半径:"+r+",改正轨道倾角:"+i
+",\n卫星在轨道上的位置("+x+","+y+"),赤纬:"+lambda
+",\n卫星在地心地固坐标系中的位置:("+x_E+","+y_E+","+z_E+")\n\n";
}
public static void main(String[] args) throws FileNotFoundException,IOException {
//打开导航电文文件 windows文件地址手动输 复制会产生空格
File file=new File("D:\\gnss\\gnss.txt");
//字符流
BufferedReader br=new BufferedReader(new FileReader(file));
String line=null;
StringBuffer navMsg=new StringBuffer();
while((line=br.readLine())!=null) {
//System.out.println(line);
navMsg.append(line+"\n");
}
//类型转换 StringBuffer->String
String navStr=navMsg.toString();
br.close();
//处理格式问题
//将字符串拆分为单个字符数组
String[] arrayChar=navStr.split("");
// for(String i:arrayChar) {
// System.out.print(i+",");
// }
ArrayList<String> al=new ArrayList<String>();
for(String element:arrayChar) {
al.add(element);
// System.out.print(element+",");
}
ListIterator<String> iter=al.listIterator();
//格式修改
while(iter.hasNext()) {
String ele=iter.next();
// System.out.print(ele);
if(ele.equals("-")) {
//后移
iter.next();
//后移
String after1=iter.next();
//判断负号后第二位是否为"."
if(after1.equals(".")) {
//前移
iter.previous();
//前移
iter.previous();
//前移
iter.previous();
//替换
iter.set(" -");
}
}
}
String[] toarray=(String[])al.toArray(new String[0]);
StringBuffer sb=new StringBuffer();
for(String element:toarray) {
sb.append(element);
// System.out.print(element);
}
String newNavStr=sb.toString();
//拆分字符串,生成数组,数组中的每一个字符串代表卫星的一个数据组
String[] arrayGroup=newNavStr.split(" ");
StringBuffer sb_output=new StringBuffer();
//头尾去掉
for(int i=2;i<arrayGroup.length-1;i++) {
String[] arrayElement=arrayGroup[i].split(" ");
// for(String j:arrayElement) {
// System.out.print(j+",");
// }
//数据龄期
String satellite=arrayElement[0].split("\n")[1];
Double IODE=Double.parseDouble(arrayElement[14]);
//Crs
Double Crs=Double.parseDouble(arrayElement[15]);
//角速度改正数
Double delta_n=Double.parseDouble(arrayElement[16]);
//参考时刻的平近点角
Double M0=Double.parseDouble(arrayElement[17]);
//Cuc
Double Cuc=Double.parseDouble(arrayElement[22]);
//偏心率
Double e=Double.parseDouble(arrayElement[23]);
//Cus
Double Cus=Double.parseDouble(arrayElement[24]);
//长半轴的开方
Double sqrt_a=Double.parseDouble(arrayElement[25]);
//toe
Double toe=Double.parseDouble(arrayElement[30]);
//Cic
Double Cic=Double.parseDouble(arrayElement[31]);
//参考时刻的升交点赤经
Double OMG0=Double.parseDouble(arrayElement[32]);
//Cis
Double Cis=Double.parseDouble(arrayElement[33]);
//参考时刻的轨道倾角
Double i0=Double.parseDouble(arrayElement[38]);
//Crc
Double Crc=Double.parseDouble(arrayElement[39]);
//参考时刻的近地点角距
Double omg=Double.parseDouble(arrayElement[40]);
//升交点赤经变化率
Double OMG0_DOT=Double.parseDouble(arrayElement[41]);
//轨道倾角变化率
Double IDOT=Double.parseDouble(arrayElement[46]);
//toc 时钟时间的周内秒
Integer year=Integer.parseInt(arrayElement[1]);
Integer month=Integer.parseInt(arrayElement[2]);
Integer date=Integer.parseInt(arrayElement[3]);
Integer hour=Integer.parseInt(arrayElement[4]);
Integer minu=Integer.parseInt(arrayElement[5]);
Integer sec=Integer.parseInt(arrayElement[6]);
LocalDateTime time=LocalDateTime.of(year, month, date, hour, minu, sec);
Date dtTime=Date.from(time.atZone(ZoneId.systemDefault()).toInstant());
Long timeSec=dtTime.getTime();
//BDT起始时间 2006年1月1日0时0分0秒 星期日
LocalDateTime time0=LocalDateTime.of(2006,1,1,0,0,0);
Date dtTime0=Date.from(time0.atZone(ZoneId.systemDefault()).toInstant());
Long time0Sec=dtTime0.getTime();
//求toc 时钟时间
Long toc=((timeSec-time0Sec)/1000)%(604800);
//钟差参数 a0,a1,a2
Double a0=Double.parseDouble(arrayElement[7]);
Double a1=Double.parseDouble(arrayElement[8]);
Double a2=Double.parseDouble(arrayElement[9]);
//设置观测时间为每分钟的第13秒
Long t0=toc+13;
LocalDateTime time_ob=LocalDateTime.of(year, month, date, hour, minu, sec+13);
Double dt=a0+a1*(t0-toc)+a2*Math.pow((t0-toc), 2);
Double t=t0-dt;
System.out.print("观测时间:"+time_ob+"\n");
System.out.print(calculateLocation(satellite,sqrt_a,delta_n,t,toe,M0,e,omg,Cuc,Cus,Crc,Crs,Cic,Cis,i0,OMG0,OMG0_DOT));
sb_output.append("观测时间:"+time_ob+"\n"+calculateLocation(satellite,sqrt_a,delta_n,t,toe,M0,e,omg,Cuc,Cus,Crc,Crs,Cic,Cis,i0,OMG0,OMG0_DOT));
}
String outputStr=sb_output.toString();
//向文件中写入计算结果
File outFile=new File("D:\\gnss\\position.txt");
FileWriter write=new FileWriter(outFile,true);
BufferedWriter rt=new BufferedWriter(write);
rt.write(outputStr);
rt.flush();
write.close();
rt.close();
}
}
适用于 RINEX 3.04格式,比较早的格式可能无法运行。(格式的原因导致我最终无法解算老师给的古早星历文件,写了个阉割版交上去了,心塞。)
放出来交流,如有错误,欢迎指正。
(感觉可以结合Cesium做个小程序玩玩。