dorvor_to_odr.m

随Doris提供的轨道数据读取程序。功能为从原始CNES ENVISAT ‘DOR_VOR ’文件中生成轨道数据记录

function dorvor_to_odr(dorvor_list)

% Function to generate Orbital Data Records (ODR) from the original
% CNES ENVISAT 'DOR_VOR' files for use with 'getorb'.
%
% Delft Institute of Earth Observation and Space Systems (DEOS), 
% Delft University of Technology, is not responsible for damage 
% of any kind caused by this script.
%
% Input:  - dorvor_list              ascii file with filenames of
%                                    DOR_VOR_* files
%
% Example input:
% DOR_VOR_AXVF-P20020906_120800_20020722_215528_20020724_002328
% DOR_VOR_AXVF-P20020906_120800_20020723_215528_20020725_002328
% DOR_VOR_AXVF-P20020906_120800_20020724_215528_20020726_002328
% DOR_VOR_AXVF-P20020906_120900_20020725_215528_20020727_002328
% DOR_VOR_AXVF-P20020906_121000_20020726_215528_20020728_002328
% DOR_VOR_AXVF-P20020906_121100_20020727_215528_20020729_002328
% DOR_VOR_AXVF-P20020906_121200_20020728_215528_20020730_002328
% DOR_VOR_AXVF-P20020912_155700_20020729_215528_20020731_002328
% ...
%
% ----------------------------------------------------------------------
% File............: dorvor_to_odr
% Version & Date..: 0.1, 21-JAN-2007
% Author..........: Freek van Leijen
%                   Delft Institute of Earth Observation and Space Systems
%                   Delft University of Technology
% ----------------------------------------------------------------------
% 

files_in = textread(dorvor_list,'%s');
date0 = datenum('01-Jan-1985');

fid1 = fopen('arclist','w');
fprintf(fid1,'%s\n','			CNES ENVISAT ORBITAL DATA RECORDS');
fprintf(fid1,'\n');
fprintf(fid1,'%s\n','This table presents some information on the CNES ENVISAT Orbital Data Records (ODR)');
fprintf(fid1,'%s\n','generated from the original ''DOR_VOR'' files by a script of the Delft Institute');
fprintf(fid1,'%s\n','of Earth Observation and Space Systems (DEOS).');
fprintf(fid1,'%s\n','DEOS, Delft University of Technology, is not responsible');
fprintf(fid1,'%s\n','for any damage of any kind caused by this script.');
fprintf(fid1,'\n');
fprintf(fid1,'%s\n','Store this file with the name "arclist" in the same directory as the ODR');
fprintf(fid1,'%s\n','files, named "ODR.???".');
fprintf(fid1,'\n');
fprintf(fid1,'%s\n','Each ODR file contains the orbital position of ENVISAT as a function of time,');
fprintf(fid1,'%s\n','computed in one orbit generation run (arc). These arcs are indicated by a');
fprintf(fid1,'%s\n','three digit number ("Arc#"). The "Arc interval" is the period over which');
fprintf(fid1,'%s\n','the orbit was computed and overlaps with the preceeding and consecutive');
fprintf(fid1,'%s\n','arc. The residuals of the tracking data (all in centimeters) are given');
fprintf(fid1,'%s\n','under "SLR", "xover", and "altim", respectively.');
fprintf(fid1,'%s\n','The length of the repeat cycle in days ("Repeat") and the distribution version');
fprintf(fid1,'%s\n','number ("Ver") are listed. The recommended begin of the precise (middle) part');
fprintf(fid1,'%s\n','of the arc ("Begin") is always over the Antarctic.');
fprintf(fid1,'\n');
fprintf(fid1,'%s\n','Arc# ------- Arc interval ------ -SLR-xover-altim  Repeat Ver  ---- Begin ----');


files_in_char = char(files_in);
months = unique(str2num(files_in_char(:,31:36)));
Nmonths = size(months,1);

for v = 1:Nmonths
  
  display(['Processing month ' num2str(v) ' of ' num2str(Nmonths) ' ...']);
  
  index = strmatch(num2str(months(v)),files_in_char(:,31:36));
  Ndays = size(index,1);
  [days,index2] = unique(str2num(files_in_char(index,31:38)));
  index = index(index2);
  Ndays = size(index,1);
  [days_sort,index_sort] = sort(days);
  
  Ndata_old = 0;
  dates = cell(100000,1);
  times = cell(100000,1);
  x = cell(100000,1);
  y = cell(100000,1);
  z = cell(100000,1);
  
  for w = 1:Ndays
    data = textread(char(files_in(index(index_sort(w)))),'%s');
    offset = strmatch('DSR_SIZE',data);
    dates_temp = data(offset+1:11:end);
    times_temp = data(offset+2:11:end);
    x_temp = data(offset+5:11:end);
    y_temp = data(offset+6:11:end);
    z_temp = data(offset+7:11:end);
    
    if w>1
      date_index = strmatch(date_end,dates_temp);
      time_index = strmatch(time_end,times_temp(date_index));
      
      Noverlap = date_index(time_index);
      if ~isempty(Noverlap)
        dates_old = dates(Ndata_old-Noverlap+1:Ndata_old);
        times_old = times(Ndata_old-Noverlap+1:Ndata_old);
        x_old = x(Ndata_old-Noverlap+1:Ndata_old);
        y_old = y(Ndata_old-Noverlap+1:Ndata_old);
        z_old = z(Ndata_old-Noverlap+1:Ndata_old);
        dates_new = dates_temp(1:Noverlap);
        times_new = times_temp(1:Noverlap);
        x_new = x_temp(1:Noverlap);
        y_new = y_temp(1:Noverlap);
        z_new = z_temp(1:Noverlap);
        
        if length(z_new)~=length(z_old)
          error('Overlap size does not match');
        end
        if isempty(strmatch(dates_new(1),dates_old(1)))
          error('The first date does not match');
        end
        if isempty(strmatch(dates_new(Noverlap),dates_old(Noverlap)))
          error('The last date does not match');
        end
        if isempty(strmatch(times_new(1),times_old(1)))
          error('The first time does not match');
        end
        if isempty(strmatch(times_new(Noverlap),times_old(Noverlap)))
          error('The last time does not match');
        end
        
        weight1 = (1:1:Noverlap)'/(Noverlap+1);
        weight2 = flipud(weight1);
        
        x_update = weight1.*str2num(char(x_new))+weight2.*str2num(char(x_old));
        y_update = weight1.*str2num(char(y_new))+weight2.*str2num(char(y_old));
        z_update = weight1.*str2num(char(z_new))+weight2.*str2num(char(z_old));
        
        x(Ndata_old-Noverlap+1:Ndata_old) = cellstr(num2str(x_update));
        y(Ndata_old-Noverlap+1:Ndata_old) = cellstr(num2str(y_update));
        z(Ndata_old-Noverlap+1:Ndata_old) = cellstr(num2str(z_update));
      end

      dates_temp(1:Noverlap) = [];
      times_temp(1:Noverlap) = [];
      x_temp(1:Noverlap) = [];
      y_temp(1:Noverlap) = [];
      z_temp(1:Noverlap) = [];
    end
    
    Ndata = size(dates_temp,1);
    
    dates(Ndata_old+1:Ndata_old+Ndata) = dates_temp;
    times(Ndata_old+1:Ndata_old+Ndata) = times_temp;
    x(Ndata_old+1:Ndata_old+Ndata) = x_temp;
    y(Ndata_old+1:Ndata_old+Ndata) = y_temp;
    z(Ndata_old+1:Ndata_old+Ndata) = z_temp;
    Ndata_old = Ndata_old+Ndata;
    
    date_end = dates_temp(end);
    time_end = times_temp(end);
    x_end = x_temp(end);
    y_end = y_temp(end);
    z_end = z_temp(end);
  end
  
  dates(Ndata_old+1:end) = [];
  times(Ndata_old+1:end) = [];
  times = char(times);
  x(Ndata_old+1:end) = [];
  y(Ndata_old+1:end) = [];
  z(Ndata_old+1:end) = [];
  
  dates_num = datenum(dates)-date0;
  seconds = dates_num*24*3600;
  hh = str2num(times(:,1:2));
  mm = str2num(times(:,4:5));
  ss = str2num(times(:,7:end));
  seconds_day = 3600*hh+60*mm+ss;
  seconds = seconds+seconds_day;
 
  xyz = [str2num(char(x)) str2num(char(y)) str2num(char(z))];
  plh = xyz2plh_trans(xyz,'GRS80_getorb'); %getorb uses GRS80
                                           %(rounded flattening)
  % getorb ODR files contain data in phi,lambda,height with
  % respect to the GRS80 ellipsoid (rounded flattening). Doris uses
  % the x,y,z output of getorb, which is equal for WGS84 and GRS80
  % (i.e., same origin and axes), so no problem there. Just have to
  % make sure ODR files are in GRS80 format.
					   
				   
  phi = (180*plh(:,1)/pi)*10^7;
  lambda = (180*plh(:,2)/pi)*10^7;
  height = plh(:,3)*1000;
  output = [seconds phi lambda height];

  fid = fopen(['ODR.' num2str(v,'%03d')],'w');
  fwrite(fid,'xODR','*char');
  fwrite(fid,'CNES_ENV','*char');
  fwrite(fid,swapbytes(int32(output(101,1))),'int32');
  fwrite(fid,swapbytes(int32(35000)),'int32');
  fwrite(fid,swapbytes(int32(v)),'int32');
  fwrite(fid,swapbytes(int32(size(output,1))),'int32');
  fwrite(fid,swapbytes(int32(-1)),'int32');
  fwrite(fid,swapbytes(int32(output')),'int32');
  fclose(fid);
  
  date_start = char(datestr(datenum(dates(1)),25));
  date_stop = char(datestr(datenum(dates(end)),25));
  date_start2 = char(datestr(datenum(dates(101)),25));
  fprintf(fid1,'%s',num2str(v,'%03d'));
  fprintf(fid1,'%s','  ');
  fprintf(fid1,'%s',[date_start(1:2) date_start(4:5) date_start(7:8)]);
  fprintf(fid1,'%s',' ');
  fprintf(fid1,'%s',[num2str(hh(1,:),'%02d') ':' num2str(mm(1,:),'%02d')]);
  fprintf(fid1,'%s',' - ');
  fprintf(fid1,'%s',[date_stop(1:2) date_stop(4:5) date_stop(7:8)]);
  fprintf(fid1,'%s',' ');
  fprintf(fid1,'%s',[num2str(hh(end,:),'%02d') ':' num2str(mm(end,:),'%02d')]);
  fprintf(fid1,'%s','                   ');
  fprintf(fid1,'%s','35.000');
  fprintf(fid1,'%s',' ');
  fprintf(fid1,'%s',' -1');
  fprintf(fid1,'%s','  ');
  fprintf(fid1,'%s',[date_start2(1:2) date_start2(4:5) date_start2(7:8)]);
  fprintf(fid1,'%s',' ');
  fprintf(fid1,'%s',[num2str(hh(101,:),'%02d') ':' num2str(mm(101,:),'%02d') ':' num2str(ss(101,:),'%02d')]);
  fprintf(fid1,'\n');
  
end

fclose(fid1);


%subfunction

function plh = xyz2plh_trans(xyz,ellips,method)
%XYZ2PLH Cartesian Coordinates to Ellipsoidal coordinates.
%        Converts cartesian coordinates X, Y and Z into ellipsoidal 
%        coordinates Phi, Lambda and h:
% 
%             plh = xyz2plh_trans(xyz,ellips)
%
%        Ellips is a text string with the name of the ellipsoid or
%        a vector with the semi-major axis a and flattening 1/f.
%        Default for ellips is 'WGS-84'.
%        This subroutine uses Bowring's method by default. The more
%        conventional iterative method can be also be used
% 
%             plh = xyz2plh_trans(xyz,ellips,1)
%
%        This method is less precise on the surface of the earth, and should
%        only be used above 10-20 km of height.

%        H. van der Marel, 07-05-95
%        (c) DEOS, Delft University of Technology


switch ellips
 case 'WGS-84'
  a = 6378137;
  f = 1/298.257223563;
  GM = 3.986005e14;
 case 'GRS80'
  a = 6378137;
  f = 1/298.257222101;
  GM = 3.986005e14;
 case 'GRS80_getorb'
  a = 6378137;
  f = 1/298.257;
  GM = 3.986005e14;
end

if nargin<3, method=0;, end

% excentricity e (squared) and semi-minor axis
e2 = 2*f - f^2;
b=(1-f)*a;

[m,n]=size(xyz);
if n~=3 & m==3, xyz=xyz';, end

r  = sqrt(xyz(:,1).^2+xyz(:,2).^2);

if method==1
% compute phi via iteration
  Np = xyz(:,3);
  for i=1:4
    phi = atan( (xyz(:,3) + e2.*Np) ./ r );
    N = a ./ sqrt(1 - e2 .* sin(phi).^2);
    Np = N .* sin(phi);
  end
else
% compute phi using B.R. Bowring's equation (default method)
  u    = atan2 ( xyz(:,3).*a , r.*b ); 
  phi  = atan2 ( xyz(:,3) + (e2/(1-e2)*b) .* sin(u).^3, r - (e2*a) .* cos(u).^3 );
  N = a ./ sqrt(1 - e2 .* sin(phi).^2);
end

plh  = [ phi                       ...
         atan2(xyz(:,2),xyz(:,1))  ...
         r ./ cos(phi) - N          ];

if n~=3 & m==3, plh=plh';, end