wuchafenxi.zip_核爆

function [Mtion,Mtropo,Mtropo1]=tion1(xy) load id75.mat; c=299792458; pi=3.1415926; for i=1:46 b(i)=id75(i); end L=46; D=fix(L/46); for j=1:D for i=5+46*(j-1):12+46*(j-1) if ( b(i)>127) b(i)= b(i)-256; else b(i)= b(i); end; end; i=5; a0(j)= b(i+46*(j-1))*2^(-28); a1(j)= b(i+1+46*(j-1))*2^(-27); a2(j)= b(i+2+46*(j-1))*2^(-24); a3(j)= b(i+3+46*(j-1))*2^(-24); b0(j)= b(i+4+46*(j-1))*2^(11); b1(j)= b(i+5+46*(j-1))*2^(14); b2(j)= b(i+6+46*(j-1))*2^(16); b3(j)= b(i+7+46*(j-1))*2^(16); %经纬度（单位弧度） fM=0.5016; lbudM=2.0222; [X,Y,Z,Xkdot,Ykdot,Zkdot,V,sv]=sat_postions(xy); D=size(X); %接收机坐标 X0=-2441605.20; Y0=5036172.91; Z0=3048605.82; zx=X-X0; zy=Y-Y0; zz=Z-Z0; %就进行一组数据进行转换就可以 H=[-sin(fM(1))*cos(lbudM(1)) -sin(fM(1))*sin(lbudM(1)) cos(fM(1));-sin(lbudM(1)) cos(lbudM(1)) 0;cos(fM(1))*cos(lbudM(1)) cos(fM(1))*sin(lbudM(1)) sin(fM(1))] ; for i=D(2):-1:1 j=j+1; A1(j)=zx(i); j=j+1; A1(j)=zy(i); j=j+1; A1(j)=zz(i); end; for i=1:24 A2(i)=A1(i+1); end; ZB=rot90(reshape(A2,3,D(2))); DXZB=ZB*H'; for i=1:D(2); DXZBX(i)=DXZB(i); DXZBY(i)=DXZB(i+D(2)); DXZBZ(i)=DXZB(i+2*D(2)); end; E=atan(DXZBZ./((DXZBX.^2+DXZBY.^2).^(0.5)))/pi; A=atan2(DXZBY,DXZBX); if(A<0) A=A+2*pi; end; fai=0.0137./(E+0.11)-0.022; fi=fM/pi+fai.*cos(A); %D=fix(l(2)/46); %for m=1:D if(abs(fi)<=0.416) fi=fi; else if(fi>0.416) fi=0.416; else fi=-0.416; end; end; lbudi=lbudM/pi+fai.*sin(A)/cos(fi*pi); fm=fi+0.064*cos((lbudi-1.617)*pi); Ag=a0+a1.*fm+a2.*(fm).^2+a3.*(fm).^3; Pg=b0+b1.*fm+b2.*(fm).^2+b3.*(fm).^3; end; for a=1:D(2) if(Ag(a)<0) Ag(a)=0; end; if(Pg(a)<72000) Pg(a)=72000; end; end; F=1+16.*((0.53-E).^3); [PGPST,SV,deltat,PRm,Numbe_measurementblocks,Carrier_Phase_Range]=Pseudoranges(xy); %计算当地时间 PGPST=11*3600-10; tion=(5*10^(-9)+Ag.*cos(2*pi./Pg.*(PGPST-14*3600))).*F; Mtion=tion*c; if(E<5/180) Mtion=200; else Mtion=Mtion; end; h=37.3; %[H]=receivers_XYZ(0); Mtropo = 2.4224.*exp(-0.13346e-3*h)./(0.026+sin(E*pi)); tropo1=0.0027*sec(E*pi).*(760+(1255/293+0.05)*15-tan(E*pi).^2)/c;%另一种算法 Mtropo1=tropo1*c;