matlab求助 讀取資料

%把每個公式里的變數都放在1×417的矩陣里，把公式里的變數代換成 用矩陣和矩陣運算，這樣就避免回圈變數的問題，直接出來417種資料，缺點就是麻煩
function Orbit = readrinexG1(File,varargin)
if size(varargin,2)==1; Prn = varargin{1};else Prn = 1:36;end

Orbit.Prn   = []; % Prn number
Orbit.Coord = []; % Cartesian coordinates,  metres
Orbit.Time  = []; % Times of coordinates,   datenums
Orbit.Clock = []; % Satellite clock errors, microsec
Orbit.Sat   = [];
Orbit       = repmat(Orbit,1,length(Prn));

pi=3.1415926;
GM=3.986005*10^14;
for k = 1:size(File,1)
        if exist(File(k,:),'file') %如果第k行存在“型別是file”的值
        x = textread(File,'%s','delimiter','\n');
        x = strvcat(x{:});
        izero = find(x(1:end,1)=='G');
        i=izero(4:end,1);
        Time= datenum(str2num(x(i,5:23)));
        Sat = str2num(x(i,2:3));
        shu=numel(i);
        %建立一個598行，每列元素是GM的矩陣
        a=ones(shu,1);
        a(a==1)=GM;
        % 讀取第1行
        year=str2num(x(i,5:8));
        month=str2num(x(i,10:11));
        day=str2num(x(i,13:14));
        hour=str2num(x(i,16:17));
        min=str2num(x(i,19:20));
        sec=str2num(x(i,22:23));
        %衛星鐘的偏差s，漂移s/s，漂移速度s/(s^2)
        af0=str2num(x(i,24:42));%√
        af1=str2num(x(i,43:61));%√
        af2=str2num(x(i,62:80));%√
        % 讀取第2行
        %星歷發布時間IODE,Crs(m),O_n(rad/s),M0(rad)  
        IODE=str2num(x(i+1,1:18));%√
        Crs=str2num(x(i+1,19:37));%√
        O_n=str2num(x(i+1,38:56)); %√
        M0=str2num(x(i+1,57:80));%√
%以上提取資料均正確
        %讀取第3行
        %衛星的Cuc(rad),軌道偏心率e,Cus(radions),Sqrt(A)(米的開方)
        Cuc=str2num(x(i+2,1:19));%√
        e=str2num(x(i+2,20:38));%√
        Cus=str2num(x(i+2,39:57));%√
        Sqrt_A=str2num(x(i+2,58:80));%√
%以上提取資料均正確 
        %讀取第4行   
        %TOE星歷的參考時刻Toe,Cic(rad),Omega0(rad),Cis(rad)
        toe=str2num(x(i+3,1:18));%√
        Cic=str2num(x(i+3,19:37));%√
        Omega0=str2num(x(i+3,38:56));%√
        Cis=str2num(x(i+3,57:80));%√
        %讀取第5行
        %參考時刻升交點赤徑i_0(rad),Crc(m),   w(rad),Omega(rad)   
        i_0=str2num(x(i+4,1:19));%√
        Crc=str2num(x(i+4,20:37));%√
        w=str2num(x(i+4,38:56));%√
        Omega=str2num(x(i+4,57:80));%√
%以上提取資料均正確  
        %讀取第6行
        %軌道傾角變化率IDOT i(rad/s),L2上的碼cflgL2，PS周數weekno，L2—P碼資料標記pflgL2  
        IDOT_i=str2num(x(i+5,1:19));%√
        % cflgL2=str2num(x(i+5,20:37));
        weekno=str2num(x(i+5,39:57));%√
        %以上提取資料均正確   
        %讀取第7行
        %衛星精度Svac(m),衛星健康狀態Svhlth,TGD(sec),IODC鐘的資料齡期  
        Svac=str2num(x(i+6,1:19));%√
        Svhlth=str2num(x(i+6,20:37));%√
        TGD=str2num(x(i+6,39:57));%√
        IDOC=str2num(x(i+6,58:80));%√ 
%以上提取資料均正確    
        %讀取第8行  這行沒用
        %電文發送時刻Transmist,擬合區間h. 
        % Transmit(i)=str2num(line(4:22)); 
       
% % % 計算衛星坐標 % % %         
%計算衛星平均角速度n
        n=(sqrt(a)./(Sqrt_A.^3))+O_n;%√
        %計算觀測瞬間衛星的平近點角M
        ut=hour+(min/60)+(sec/3600);%先將民用日的時分秒化為實數時%√
%         if month<=2
%         year=year-1;
%         month=month+12;
%         else
%         year=year;
%         month=month;
%         end
        % if month(i)<=2
        % y(i)=year(i)-1;
        % m(i)=month(i)+12;
        % else
        % y(i)=year(i);
        % m(i)=month(i);
        % end
        jd=fix(365.25*year)+fix(30.6001*(month+1))+day+(ut/24)+1720981.5;%將民用日的年月日轉換成儒略日%√
        % weekgps(i)=int((jd(i)-2444244.5)/7);    
        weeksec=(jd-2444244.5-weekno*7)*24*3600;   %√
        tk=weeksec-toe;%？？？？？
        deltat=af0+af1.*tk+af2.*tk.^2;%？？？？
        tk=tk+deltat;  %？？？？？
      
        sdf=find(tk>302400);%√
        tk(sdf,:)=tk(sdf,:)-604800;%√
        sdf1=find(tk<-302400);  %√
        tk(sdf1,:)=tk(sdf1,:)+604800;%√
        M=M0+n.*tk;%√
        % % 計算偏近點角
        E_temp=M;%√
        d_E=1;%√
while(abs(d_E)<1e-15) %？？？？？
        E_temp_next=M+e.*sin(E_temp);%？？？？？
        d_E=E_temp_next-E_temp;%？？？？？
        E_temp=E_temp_next;%？？？？？
end
        E=E_temp;
        V=atan2(sqrt(1-e.*e).*sin(E),(cos(E)-e));%√

        %計算升交距u
        u=V+w;
        %計算攝動改正項
        SigmaU=Cuc.*cos(2*u)+Cus.*sin(2*u);
        SigmaR=Crc.*cos(2*u)+Crs.*sin(2*u);
        SigmaI=Cic.*cos(2*u)+Cis.*sin(2*u);
        %對u,r,i,進行攝動改正
        u=u+SigmaU;
        a=(Sqrt_A).^2;
        r=a.*(1-e.*cos(E))+SigmaR;
        I=i_0+SigmaI+IDOT_i.*tk;
        %計算衛星在軌道坐標系中的位置
        x0=r.*cos(u);
        y0=r.*sin(u);
        % xy(1,ii)=x0;
        % xy(2,ii)=y0;
        %計算觀測瞬間升交點的經度L
        we1=7.29211567*10^(-5);%地球自轉角速度
        we=ones(shu,1);
        we(we==1)=we1;
        L=Omega0+(Omega-we).*tk-we.*toe;
        %計算衛星在瞬時地球坐標系中的位置
        %XYZ=struct(i'X',NaN,'Y',NaN,'Z',NaN);
        X=x0.*cos(L)-y0.*cos(I).*sin(L);
        Y=x0.*sin(L)+y0.*cos(I).*cos(L);
        Z=y0.*sin(I);
        Pos=[X Y Z]; 
        
        for i = 1:length(Prn)
        Orbit(i).Prn   = Prn(i);          
        Orbit(i).Time  = [Orbit(i).Time  ;Time(find(Sat==Prn(i)),:)];
        Orbit(i).Coord = [Orbit(i).Coord ;Pos(find(Sat==Prn(i)),:)];
%         Orbit(i).Clock = [Orbit(i).Clock ;Clk(find(Sat==Prn(i)),:)]; 
        Orbit(i).Sat   = [Orbit(i).Sat ;Sat(find(Sat==Prn(i)),:)]; 
        end
         else fprintf('Warning from ning: %s not found\n',File(k,:)); end
end
%因為有時間重復的時候，所以后面Coord資料也有重復的，后面程式自動洗掉重復的數
% 所以原本20個資料，結果為11個，剩下九個是重復的
for i = 1:length(Orbit)
      if ~isempty(Orbit(i).Coord)
       
      [t,j] = unique(round(Orbit(i).Time*24*60));
      %  [t,j] = unique(round(Orbit(i).Time));
      Orbit(i).Time  = Orbit(i).Time(j);
      Orbit(i).Coord = Orbit(i).Coord(j,:);
%       Orbit(i).Clock = Orbit(i).Clock(j);
   end
end

uj5u.com熱心網友回復：

在斷點處   讀取不到   20-38字符   讀取不到  小 e  矩陣

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/202807.html

標籤：其他開發語言

上一篇：Python機器學習演算法：線性回歸

下一篇：使用LIME解釋黑盒ML模型