根据北斗广播星历计算卫星位置的小程序 Java版本 RINEX VERSION3.04

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做个小程序玩玩。

  • 2
    点赞
  • 26
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值