EASY8.M
上传用户:sfyaiting
上传日期:2009-10-25
资源大小:320k
文件大小:3k
源码类别:

GPS编程

开发平台:

Matlab

  1. % EASY8    Test for cycle slip and repair of receiver clock offset. 
  2. %          We use the ionosphere free linear combination.
  3. %          Formulas are described in Strang and Borre, page 491
  4. %Kai Borre 27-07-2002
  5. %Copyright (c) by Kai Borre
  6. %$Revision: 1.0 $  $Date: 2002/07/27  $
  7. v_light = 299792458;    % vacuum speed of light m/s
  8. f0 = 10.23*10^6;
  9. f1 = 154*f0;
  10. f2 = 120*f0;
  11. lambda1 = v_light/f1;
  12. lambda2 = v_light/f2;
  13. alpha1 = f1^2/(f1^2-f2^2);
  14. alpha2 = 1-alpha1;
  15. ofile = 'site247j.01o';
  16. %ofile = 'kofi2.01o';
  17. fid = fopen(ofile,'rt');
  18. [Obs_types, ant_delta, ifound_types, eof1] = anheader(ofile);
  19. if ((ifound_types == 0) | (eof1 == 1))
  20.     error('Basic information is missing in RINEX file')
  21. end;
  22. NoObs_types = size(Obs_types,2)/2;
  23. j = fobs_typ(Obs_types,'P1');
  24. k = fobs_typ(Obs_types,'P2');
  25. l = fobs_typ(Obs_types,'L1');
  26. m = fobs_typ(Obs_types,'L2');
  27. cols = [j k l m];
  28. fid = fopen(ofile,'rt');
  29. for i = 1:20
  30.     [tr_RAW, dt, sv, eof2] = fepoch_0(fid); 
  31.     NoSv = size(sv,1);
  32.     obs = grabdata(fid, NoSv, NoObs_types);  
  33.     Obs = obs(:,cols);
  34.     for j = 1:NoSv
  35.         if i == 1
  36.             Obs0 = Obs;
  37.             P(:,1) = alpha1*Obs0(:,1)+alpha2*Obs0(:,2);
  38.             Phi(:,1) = alpha1*lambda1*Obs0(:,3)+alpha2*lambda2*Obs0(:,4);
  39.         else
  40.             P(j,i) = alpha1*Obs(j,1)+alpha2*Obs(j,2);
  41.             deltaP(j,i) = P(j,i)-P(j,1);
  42.             Phi(j,i) = alpha1*lambda1*Obs(j,3)+alpha2*lambda2*Obs(j,4);
  43.             deltaPhi(j,i) = Phi(j,i)-Phi(j,1);
  44.         end
  45.     end
  46. end
  47. fclose all;
  48. DP = diff(deltaP,1,2);
  49. % Repair of clock reset of 1ms ~ 299 km; affects only pseudoranges
  50. i1 = find(abs(DP(1,:)) > 280000);
  51. for j = i1
  52.    if DP(:,j) < 0
  53.        DP(:,j) = DP(:,j)+299792.458; 
  54.    else 
  55.        DP(:,j) = DP(:,j)-299792.458; 
  56.    end
  57. end
  58. DPhi = diff(deltaPhi,1,2);
  59. DDP = diff(DP,1,2);
  60. DDPhi = diff(DPhi,1,2);
  61. figure(1);
  62. subplot(2,2,1), plot(DP')
  63. set(get(gca,'ylabel'),'String',({'Pseudoranges Differenced','over Time [m]'}),...
  64.     'fontsize',12,'color','r','verticalalignment','bottom','rotation',90)   
  65. subplot(2,2,2), plot(DPhi')
  66. set(get(gca,'ylabel'),'String',({'Phases Differenced','over Time [m]'}),...
  67.     'fontsize',12,'color','r','verticalalignment','bottom','rotation',90)   
  68. subplot(2,2,3), plot(DDP')
  69. set(get(gca,'ylabel'),'String',({'Double Differenced','Pseudoranges [m]'}),...
  70.     'fontsize',12,'color','r','verticalalignment','bottom','rotation',90)   
  71. subplot(2,2,4), plot(DDPhi')
  72. set(get(gca,'ylabel'),'String',({'Double Differenced','Phases [m]'}),...
  73.     'fontsize',12,'color','r','verticalalignment','bottom','rotation',90)   
  74. %ax = axes('units','normal','position',[.075 .075 .9 .9],'visible','off');
  75. %set(get(ax,'title'),'visible','on')
  76. %title('Check of Cycle Slips','fontsize',16,'fontweight','bold')
  77. %h = get(ax,'title');
  78. print -deps easy8
  79. %%%%%%%%% end easy8.m %%%%%%%%%