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');  endfclose(fid1);%subfunctionfunction 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 Technologyswitch 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;endif nargin<3, method=0;, end% excentricity e (squared) and semi-minor axise2 = 2*f - f^2;b=(1-f)*a;[m,n]=size(xyz);if n~=3 & m==3, xyz=xyz';, endr  = 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);  endelse% 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);endplh  = [ phi                       ...         atan2(xyz(:,2),xyz(:,1))  ...         r ./ cos(phi) - N          ];if n~=3 & m==3, plh=plh';, end