function []=ps_export_csv(ts_type,ifg_list,filename)

%function []=ps_export_csv(velocity_type,ts_type,filename)
% PS_EXPORT_CSV export velocity and time series data to csv file. The mean
%       velocities are in mm/year, the unwrapped time series are cumulative
%       displacements in mm. In addition the DEM error and the coherence
%       estimates are exported in a separate column.
%    PS_EXPORT_CSV(VELOCITY_TYPE,TS_TYPE,FILENAME)
%
%    valid VELOCITY_TYPE's are:
%      
%       'v', the mean velocity from single master interferograms
%
%       'vsb', the mean velocity from sbas interferograms
%       
%        'v-dmos', the mean velocity from single master interferograms
%           corrected by the DEM error estimated in the time domain,
%           orbital errors as well as spatially correlated noise
%       
%       'v-dsbmos', the mean velocity from single master interferograms
%           corrected by the DEM error in the interferogram domain
%
%
%   valid TS_TYPES are: 
%
%       'u-dmos' unwrapped phase corrected by the DEM error estimated in
%           the time domain, master atmosphere, orbital errors and slave
%           atmosphere.
%   
%       'u-dsbmos' unwrapped phase corrected by the DEM error estimated in
%           interferogram domain.
%
%       'd' cumulative displacement that would result from the DEM error
%           estimated in the time domain
%
%       'dsb' cumulative displacement that would result from the DEM error
%           estimated in the interferogram domain

stdargin = nargin ; 

if stdargin<3
    help ps_export_csv
    error('not enough input args')
end

%for testing puposes:
% velocity_type='v-dmos';
%ts_type='u-dsbmos';
%ifg_list=[];
% 
%filename=name_d;
%% Load the names of the matlab files
load psver
bpname=['./bp',num2str(psver)];
psname=['./ps',num2str(psver)];
pmname=['./pm',num2str(psver)];
rcname=['./rc',num2str(psver),'.mat'];
phname=['./ph',num2str(psver),'.mat'];
phuwname=['./phuw',num2str(psver)];
phuwsbname=['./phuw_sb',num2str(psver)];
phuwsbresname=['./phuw_sb_res',num2str(psver)];
scnname=['./scn',num2str(psver)];
ifgstdname=['./ifgstd',num2str(psver)];
apsbandsname = ['./tca_bands' num2str(psver) '.mat'];
apsbandssbname = ['./tca_bands_sb' num2str(psver) '.mat'];
apsname=['./tca',num2str(psver)];
apssbname=['./tca_sb',num2str(psver)];
hgtname=['./hgt',num2str(psver)];
sclaname=['./scla',num2str(psver)];
sclasbname=['./scla_sb',num2str(psver)];
sclasmoothname=['./scla_smooth',num2str(psver)];
sclasbsmoothname=['./scla_smooth_sb',num2str(psver)];
meanvname=['./mv',num2str(psver)];
filtername=['./filt_sb',num2str(psver)];


iononame = ['./ica' num2str(psver) '.mat'];
ionosbname = ['./ica_sb' num2str(psver) '.mat'];
tidename=['./tide',num2str(psver)];
tidesbname=['./tide_sb',num2str(psver)];


%% Get various parameters
drop_ifg_index=getparm('drop_ifg_index');
small_baseline_flag=getparm('small_baseline_flag');
scla_deramp = getparm('scla_deramp');
lambda=getparm('lambda');

%% Get time info

ps=load(psname);
day=ps.day;
master_day=ps.master_day;
xy=ps.xy;
lonlat=ps.lonlat;
n_ps=ps.n_ps;
n_ifg=ps.n_ifg;
master_ix=sum(day<master_day)+1;
ref_ps=0;    
%% Start switching through the velcocity types

switch(ts_type)
    case {'u-dmos'}
        ref_ps=ps_setref;
        disp('we correct the time series with single master dem error');
        uw=load(phuwname);
        scn=load(scnname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        scla=load(sclaname);
        ts=uw.ph_uw - scn.ph_scn_slave - repmat(scla.C_ps_uw,1,size(uw.ph_uw,2)) - scla.ph_scla - scla.ph_ramp;
        ts(:,ps.master_ix)=0;
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        else
            unwrap_ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);
        end
        clear uw scn scla
        
        
 
        
        
        case {'u-mos'}
        ref_ps=ps_setref;
        disp('we correct the time series with single master dem error');
        uw=load(phuwname);
        scn=load(scnname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        scla=load(sclaname);
        ts=uw.ph_uw - scn.ph_scn_slave - repmat(scla.C_ps_uw,1,size(uw.ph_uw,2)) - scla.ph_ramp;
        ts(:,ps.master_ix)=0;
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        else
            unwrap_ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);
        end
        clear uw scn scla
        
     case {'u-dsbof'}
        ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
        filt=load(filtername);
        ph_filt=filt.ph_uw_high;
        G=G(:,unwrap_ifg_index);
        bperp_some=[G\double(bp.bperp_mat')]'; 
        bperp_mat(:,unwrap_ifg_index)=bperp_some;
        ph_scla=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;
         % Get the DEM error from the scaler
         
        ts=uw.ph_uw - ph_scla - scla.ph_ramp-ph_filt;
        clear uw scn scla
        ts(:,ps.master_ix)=0;
        
        %% Correct time series by DEM error
    case {'u'}
        ref_ps=ps_setref;
        disp('we do not correct the time series');
        uw=load(phuwname);
         ts=uw.ph_uw;
         
    case {'u-o'}
        ref_ps=ps_setref;
        disp('we do not correct the time series');
        
        uw=load(phuwname);
        
         if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
       scla=load(sclaname);
         ts=uw.ph_uw- scla.ph_ramp;
    case {'u-dsbmos'}
        ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        scn=load(scnname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
          G=G(:,unwrap_ifg_index);
          bperp_some=[G\double(bp.bperp_mat')]'; 
          bperp_mat(:,unwrap_ifg_index)=bperp_some;
          ph_scla=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;
         % Get the DEM error from the scaler
        ts=uw.ph_uw - scn.ph_scn_slave - repmat(scla.C_ps_uw,1,size(uw.ph_uw,2)) - ph_scla - scla.ph_ramp;
        clear uw scn scla
        ts(:,ps.master_ix)=0;
        
        
   case {'u-dsbmo'}
        ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
          G=G(:,unwrap_ifg_index);
          bperp_some=[G\double(bp.bperp_mat')]'; 
          bperp_mat(:,unwrap_ifg_index)=bperp_some;
          ph_scla=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;
         % Get the DEM error from the scaler
        ts=uw.ph_uw - repmat(scla.C_ps_uw,1,size(uw.ph_uw,2)) - ph_scla - scla.ph_ramp;
        clear uw scn scla
        ts(:,ps.master_ix)=0;
        
    case {'u-dsbo'}
        ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
          G=G(:,unwrap_ifg_index);
          bperp_some=[G\double(bp.bperp_mat')]'; 
          bperp_mat(:,unwrap_ifg_index)=bperp_some;
          ph_scla=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;
         % Get the DEM error from the scaler
        ts=uw.ph_uw - ph_scla - scla.ph_ramp;
        clear uw scn scla
        ts(:,ps.master_ix)=0;
        
    case {'bperp'}
        scla_drop_index=getparm('scla_drop_index');
        ref_ps=ps_setref;
        uw=load(phuwname);
        strcmpi(small_baseline_flag,'y')
        bperp_mat=zeros(ps.n_ps,ps.n_image,'single');
        G=zeros(ps.n_ifg,ps.n_image);
        for i=1:ps.n_ifg
             G(i,ps.ifgday_ix(i,1))=-1;
             G(i,ps.ifgday_ix(i,2))=1;
        end
        if isfield(uw,'unwrap_ifg_index_sm')
            unwrap_ifg_index=setdiff(uw.unwrap_ifg_index_sm,scla_drop_index);
        else
        end
        G=G(:,unwrap_ifg_index);
        ts=[G\double(bp.bperp_mat')]';

    case {'dsbts'}
        ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
          G=G(:,unwrap_ifg_index);
          bperp_some=[G\double(bp.bperp_mat')]'; 
          bperp_mat(:,unwrap_ifg_index)=bperp_some;
          ts=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;

    case{'dts'}
        scla=load(sclaname);
        ts=scla.ph_scla + scla.ph_ramp;
        
    
    case{'s'}
        scn=load(scnname);
        ts= scn.ph_scn_slave;
    
    case{'uds-dsbmos'}
        disp('We project everything downslope');
         ref_ps=ps_setref;
        disp('we correct the time series with sbas dem error');
        uw=load(phuwname);
        scn=load(scnname);
        morph=load('morph2');
                if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        sclasb=load(sclasbname);
        ph_scla_sb=sclasb.ph_scla;
        scla=load(sclaname);
        if strcmp('n',scla_deramp)
            disp('Warning: scla_deramp flag set to n. Set to y and rerun Step 7 before using the -o plot command.')
            return;
        end
        if strcmpi(small_baseline_flag,'y')
            unwrap_ifg_index=uw.unwrap_ifg_index_sm;
        end
        % Invert the scla
        G=zeros(ps.n_ifg,ps.n_image); 
        for i=1:ps.n_ifg
          G(i,ps.ifgday_ix(i,1))=-1;
          G(i,ps.ifgday_ix(i,2))=1;
        end
        bp=load(bpname);
          G=G(:,unwrap_ifg_index);
          bperp_some=[G\double(bp.bperp_mat')]'; 
          bperp_mat(:,unwrap_ifg_index)=bperp_some;
          ph_scla=repmat(sclasb.K_ps_uw,1,size(bperp_mat,2)).*bperp_mat;
         % Get the DEM error from the scaler
          ts_norm=uw.ph_uw - scn.ph_scn_slave - repmat(scla.C_ps_uw,1,size(uw.ph_uw,2)) - ph_scla - scla.ph_ramp;
          scaler=repmat(morph.scaler,1,size(ts_norm,2));
          ts=ts_norm.*scaler;
end

% subset the right interferograms
interferogram_domain='y';
if strcmpi(interferogram_domain,'y') %SBAS
    phuwres=load(phuwsbresname,'sm_cov');
    if isfield(phuwres,'sm_cov');
        sm_cov=phuwres.sm_cov(unwrap_ifg_index,unwrap_ifg_index);
    else
        sm_cov=eye(length(unwrap_ifg_index));
    end
else
    if ~exist([ifgstdname,'.mat',],'file')
        sm_cov=eye(length(unwrap_ifg_index));
    else
        ifgstd=load(ifgstdname);
        if isfield(ifgstd,'ifg_std');
            ifgvar=(ifgstd.ifg_std*pi/181).^2;
            sm_cov=diag(ifgvar(unwrap_ifg_index));
        else
            sm_cov=eye(length(unwrap_ifg_index));
        end
    end
end

 if ~isempty(ifg_list)
      unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
      sm_cov=sm_cov(unwrap_ifg_index,unwrap_ifg_index);
 end

ts=ts(:,unwrap_ifg_index);

day=ps.day(unwrap_ifg_index,:);
ts=ts-repmat(ts(:,1),1,size(ts,2));
ts=ts-repmat(nanmean(ts(ref_ps,:),1),n_ps,1);

days_referenced=day-min(day);
years_referenced=days_referenced/365.25;

%% Smooth for the first time
ph_ds_mm=ts/4/pi*(lambda)*1000;

for i=1:size(ph_ds_mm,1)
        sm=ph_ds_mm(i,:);
        m=(years_referenced\sm');
        mean_v(i)= m(1);        
end

ph_ds_mm=ph_ds_mm-repmat(nanmean(ph_ds_mm(ref_ps,:),1),n_ps,1);
ph_ds_mm=ph_ds_mm-repmat(ph_ds_mm(:,1),1,size(ph_ds_mm,2));

if strcmp(ts_type,'uds-dsbmos')
    ph_ds_mm=abs(ph_ds_mm);
    mean_v=abs(mean_v);
end


if ~strcmpi(ts_type,'bperp')
    disp('We export time series');
else
    disp('We just want the baselines');
    ph_ds_mm=ts;
end

%% Get DEM errors and coherence values
scla=load(sclaname);
if strcmpi(small_baseline_flag,'y')
    sclasb=load(sclasbname);
    demifg_m=K2q(sclasb.K_ps_uw);
    demts_m = K2q(scla.K_ps_uw);
else
    
demts_m = K2q(scla.K_ps_uw);
demifg_m=zeros(ps.n_ps,1);
end
pm=load(pmname);

coherence=pm.coh_ps;
data = horzcat(lonlat,mean_v',demts_m,demifg_m,coherence,ph_ds_mm);
 % Prepare the header
 %
formatOut = 'yyyymmdd';
dates = datestr(day,formatOut);
prefix = {'d'};
header = strcat(prefix,dates);
headerhorz= header.';
rest = {'x','y','v','dem_ts','dem_ifg','coherence'};
header = horzcat(rest,headerhorz);

disp('we write to the file')
disp(filename);
   
fid=fopen(filename,'wt');
[rows,cols]=size(header);

% Run loop to concatenate
for i=1:rows
      fprintf(fid,'%s,',header{i,1:end-1});
      fprintf(fid,'%s\n',header{i,end});
end
fclose(fid);

%% Write to file
dlmwrite(filename,data, 'precision', 10,'-append')   

end
