ps_output.m

来自「StaMps最新测试版」· M 代码 · 共 121 行

function []=ps_output()
%PS_OUTPUT write various output files 
%
%   Andy Hooper, June 2006

load psver
psname=['ps',num2str(psver)];
rcname=['rc',num2str(psver)];
phuwname=['phuw',num2str(psver)];
sclaname=['scla',num2str(psver)];
scnname=['scn',num2str(psver)];
meanvname=['mean_v'];

ps=load(psname);
phuw=load(phuwname);
%scn=load(scnname);
rc=load(rcname);

ijname=['ps_ij.txt'];
ij=ps.ij(:,2:3);
save(ijname,'ij','-ASCII');

llname=['ps_ll.txt'];
lonlat=ps.lonlat;
save(llname,'lonlat','-ASCII');

datename=['date.txt'];
date_out=str2num(datestr(ps.day,'yyyymmdd'));
save(datename,'date_out','-ASCII','-DOUBLE');

master_ix=sum(ps.master_day>ps.day)+1;

%ph_filt=phuw.ph_uw-scn.ph_scn;
ref_ps=ps_setref;
ph_uw=phuw.ph_uw-repmat(mean(phuw.ph_uw(ref_ps,:)),ps.n_ps,1);
%ph_filt=ph_filt-repmat(mean(ph_filt(ref_ps,:)),ps.n_ps,1);
ph_w=angle(rc.ph_rc.*repmat(conj(sum(rc.ph_rc(ref_ps,:))),ps.n_ps,1));
ph_w(:,master_ix)=0;


fid=fopen('ph_w.flt','w');
fwrite(fid,ph_w','float');
fclose(fid);

fid=fopen('ph_uw.flt','w');
fwrite(fid,ph_uw','float');
fclose(fid);

%fid=fopen('ph_filt.flt','w');
%fwrite(fid,ph_filt','float');
%fclose(fid);

scla=load(sclaname);
ph_uw=phuw.ph_uw - scla.ph_scla - repmat(scla.C_ps_uw,1,ps.n_image);

%%% this is only approximate
K_ps_uw=scla.K_ps_uw-mean(scla.K_ps_uw);
dem_error=double(K2q(K_ps_uw,ps.ij(:,3)));
dem_sort=sort(dem_error);
min_dem=dem_sort(ceil(length(dem_sort)*0.001));
max_dem=dem_sort(floor(length(dem_sort)*0.999));
dem_error(dem_error<min_dem)=min_dem;
dem_error(dem_error>max_dem)=max_dem;
dem_error=[ps.lonlat,dem_error];
%dem_error=[ps.lonlat,double(K_ps_uw)];

save('dem_error.xy','dem_error','-ASCII');

%%%

clear scla phuw
ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),ps.n_ps,1);

if ps.n_ifg==size(ph_uw,2)
    unwrap_ifg_index=getparm('unwrap_ifg_index')
    if strcmp(unwrap_ifg_index,'all')
        unwrap_ifg_index=[1:ps.n_ifg];
    end	
    ph_uw_sub=ph_uw(:,unwrap_ifg_index);
    day=ps.day(unwrap_ifg_index);
else
    ph_uw_sub=ph_uw;
    day=ps.day;
end
% Each ifg has master APS - slave APS, including master 
% (where slave APS = master APS) so OK to include master in inversion
%G=[ones(size(day)),day-ps.master_day]; 
%G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 

%m=G\ph_uw_sub';
meanv=load(meanvname);
lambda=getparm('lambda');
mean_v=-meanv.m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
%mean_v=-meanv.m(2,:)'*365.25/4/pi*lambda*100; % m(1,:) is master APS + mean deviation from model
v_sort=sort(mean_v);
min_v=v_sort(ceil(length(v_sort)*0.001))
max_v=v_sort(floor(length(v_sort)*0.999))
mean_v(mean_v<min_v)=min_v;
mean_v(mean_v>max_v)=max_v;
%mean_v=mean_v-min_v;

%ph_all=m(3,:)'/4/pi*lambda*1000;
%ph_all=m(1,:)'/4/pi*lambda*1000;

mean_v_name=['ps_mean_v.xy'];
mean_v=[ps.lonlat,double(mean_v)];
save(mean_v_name,'mean_v','-ASCII');

for i=1:ps.n_image
    ph=ph_uw(:,i);
    ph_sort=sort(ph);
    min_ph=ph_sort(ceil(length(ph_sort)*0.001));
    max_ph=ph_sort(floor(length(ph_sort)*0.999));
    ph(ph<min_ph)=min_ph;
    ph(ph>max_ph)=max_ph;
    ph=-ph*lambda*100/4/pi;
    ph=[ps.lonlat,double(ph)];
    save(['ps_u-dm.',num2str(i),'.xy'],'ph','-ASCII');
end