B01_ValidateTwobyTwo.m

function []=B01_ValidateTwobyTwo(fuentes,npuntos,archivo)
fuente1=fuentes{1,1};
fuente2=fuentes{2,1};

% reading information
eval(['datos_' fuente1 '=readmatrix(archivo,"Sheet","' fuente1 '", "Range","A4");'])
eval(['datos_' fuente2 '=readmatrix(archivo,"Sheet","' fuente2 '", "Range","A4");'])

%changing date
for ifuentes=1:size(fuentes,1)
    eval(['datosanalizar=datos_' fuentes{ifuentes} ';']);
    if min(min(datosanalizar(:,1:2:npuntos*2)))<10000
        anno=fix(datosanalizar(:,1:2:npuntos*2));
        bisiestoSI=leapyear(anno);
        bisiestoNO=~leapyear(anno);
        datosanalizar(:,1:2:npuntos*2)=datenum(anno,1,1)+366*bisiestoSI.*(datosanalizar(:,1:2:npuntos*2)-anno)+365*bisiestoNO.*(datosanalizar(:,1:2:npuntos*2)-anno);
        eval(['datos_' fuentes{ifuentes} '=datosanalizar;']);
    else
        datosanalizar(:,1:2:npuntos*2)=x2mdate(datosanalizar(:,1:2:npuntos*2));
        eval(['datos_' fuentes{ifuentes} '=datosanalizar;']);
    end
end

%% Detection of more than one observations in any source for the same date
for ifuentes=1:size(fuentes,1)
    eval(['datosanalizar=datos_' fuentes{ifuentes} ';']);
    datosequal=(datosanalizar(:,1:2:npuntos*2));
    tempequal= datosequal (2:end,:)- datosequal (1:end-1,:);
    [filas,col] = find (tempequal==0);
    
    if sum(filas)>1
        answer = questdlg(['Duplicate dates have been detected in source ' fuentes{ifuentes} ...
            ', is it valid for the dataset?'], 'DUPLICATE DATE', 'Yes','No','Yes');
        if strcmp (answer, 'Yes')
            for ifilas=1:length(filas)
                %if there are duplicated values, the mean is calculated and this value is assigned to the first date with data
                %media
                tempvalores = (datosanalizar(filas(ifilas),2* col(ifilas))+datosanalizar(filas(ifilas)+1, 2*col(ifilas)))/2; %mean calculation
                datosanalizar (filas(ifilas), 2*col(ifilas)) = tempvalores; %assigns to the first matching date
            end
            clear tempvalores;
            
            %remove the second redundant date and its value
            for inpuntos=1:npuntos  %
                temp1=find(col==inpuntos);
                tempdatosanalizar=datosanalizar(:,2*inpuntos-1:2*inpuntos);
                tempdatosanalizar(filas(temp1)+1,:)=[];
                tempdatosanalizar=[tempdatosanalizar;nan(length(temp1),2)];
                if inpuntos==1
                    datosanalizar2=tempdatosanalizar;
                else
                    datosanalizar2=[datosanalizar2 tempdatosanalizar];
                end
                clear temp1 tempdatosanalizar
            end
            eval(['datos_' fuentes{ifuentes} '=datosanalizar2;']);
            clear datosanalizar2 datosanalizar
        else
            error('There exist at least one duplicate date');
        end
    end
end



eval(['datoscomp1=datos_' fuente1 ';']);
eval(['datoscomp2=datos_' fuente2 ';']);



%%
angulo1=readmatrix (archivo, 'sheet',fuente1,'range','B1:B1');
angulo1=angulo1(1); 
angulo2= readmatrix (archivo, 'sheet',fuente2,'range','B1:B1');
angulo2=angulo2(1); 
col=2:2:2*npuntos;
datoscomp1(:,col)=datoscomp1(:,col)./cosd(angulo1);
datoscomp2(:,col)=datoscomp2(:,col)./cosd(angulo2);


%% Validado 
veloParcial=nan(npuntos,2);  %only two columns for the only two sources of information
estadisPun=nan(19,npuntos);
for ipuntos=1:npuntos
    col=1+(ipuntos-1)*2;
    tempdatos1=datoscomp1(:,col:col+1);  %save in "tempdatos" observed data for each observation point
    tempdatos2=datoscomp2(:,col:col+1);
    [tempf,~]=find(isnan(tempdatos1));   %removing all nan
    tempdatos1(tempf,:)=[];
    [tempf,~]=find(isnan(tempdatos2));   %removing all nan
    tempdatos2(tempf,:)=[];
    
    if ~isempty(tempdatos1) && ~isempty(tempdatos2)  %avoid comparing points without data
        fechaini=max(tempdatos1(1,1),tempdatos2(1,1));
        fechafin=min(tempdatos1(end,1),tempdatos2(end,1));
        
        if fechaini<fechafin %Check for overlap period for boths technologies
            %detecta cual comienza primero
            if tempdatos1(1,1)<tempdatos2(1,1)
                iniciaprimero=1;
            else
                iniciaprimero=2;
            end
            
            % removing the initial period without paired data from both sources
            [tempf,~]=find(tempdatos1(:,1)<fechaini);
            if ~isempty(tempf)
                tempf(end)=[];   %keep the last to delete 
                tempdatos1(tempf,:)=[];
            end
            [tempf,~]=find(tempdatos2(:,1)<fechaini);
            if ~isempty(tempf)
                tempf(end)=[];  %keep the last to delete 
                tempdatos2(tempf,:)=[];
            end
            %removing the final period without paired data from the both sources
            [tempf,~]=find(tempdatos1(:,1)>fechafin);
            if ~isempty(tempf)
                tempf(1)=[];    %keep the first to delete 
                tempdatos1(tempf,:)=[];
            end
            [tempf,~]=find(tempdatos2(:,1)>fechafin);
            if ~isempty(tempf)  %keep the first to delete 
                tempf(1)=[];
                tempdatos2(tempf,:)=[];
            end
            % 
            % From this point, "tempdatos" is only for the common period
            
            
            % Average velocity computation for the common period
            a=fitlm(tempdatos1(:,1),tempdatos1(:,2));
            veloParcial(ipuntos,1)=a.Coefficients{2,1}*365.25; %in cm/year
            a=fitlm(tempdatos2(:,1),tempdatos2(:,2));
            veloParcial(ipuntos,2)=a.Coefficients{2,1}*365.25; %in cm/year
            
            
            %% Absolute deformation validation
            %Detecting duplicate dates
            tempdatos1=sortrows(tempdatos1,1);      %sort from oldest to newest  
            tempdatos2=sortrows(tempdatos2,1);
            
            if size(tempdatos1,1)>size(tempdatos2,1)
                seriedensa=1;           %The most dense series used for interpolation
                tempdatos1prima=interp1(tempdatos1(:,1),tempdatos1(:,2),tempdatos2(:,1));
            else
                seriedensa=2;
                tempdatos2prima=interp1(tempdatos2(:,1),tempdatos2(:,2),tempdatos1(:,1));
            end
            % Next is for correcting boundary erros producing by 
            % missing data in the lineal estimation.
            % 
            if seriedensa==1
                tempdatos1mt=tempdatos1prima;
                tempdatos2mt=tempdatos2(:,2);
                tempfilaseliminar=find(isnan(tempdatos1prima));
            else
                tempdatos2mt=tempdatos2prima;
                tempdatos1mt=tempdatos1(:,2);
                tempfilaseliminar=find(isnan(tempdatos2prima));
            end
            %Filtered time series with common dates
            tempdatos2mt(tempfilaseliminar)=[];
            tempdatos1mt(tempfilaseliminar)=[];
            
            %Translation of the series for statistic computations
            if iniciaprimero==1
                tempdatos2mt=tempdatos2mt-(tempdatos2mt(1)-tempdatos1mt(1));
            elseif iniciaprimero==2
                tempdatos1mt=tempdatos1mt-(tempdatos1mt(1)-tempdatos2mt(1));
            end
            
            %save in a matrix of 3 columns: time, obs1, obs2
            if seriedensa==1
                tempdatosAmbos=tempdatos2(:,1);
            else
                tempdatosAmbos=tempdatos1(:,1);
            end
            tempdatosAmbos(tempfilaseliminar)=[];
            tempdatosAmbos=[tempdatosAmbos tempdatos1mt tempdatos2mt];
            clear tempfilaseliminar tempdatos1mt tempdatos2mt;
            
            %Calculation of absolute deformation statistics for ipuntos
            estadisPun(1,ipuntos)=max(tempdatosAmbos(:,2))-min(tempdatosAmbos(:,3));  %observed ranges
            estadisPun(2,ipuntos)=abs(mean (tempdatosAmbos(:,2)));   %absolute mean of observed points
            estadisPun(3,ipuntos)=size(tempdatosAmbos,1);           %nº of observations
            if size(tempdatosAmbos,1)>10
                estadisPun(4,ipuntos)= prctile(tempdatosAmbos(:,2),90); %Percentil 90
                estadisPun(5,ipuntos)= prctile(tempdatosAmbos(:,2),10); %Percentil 10
            else
                estadisPun(4,ipuntos)= max(tempdatosAmbos(:,2)); %Percentil 90
                estadisPun(5,ipuntos)= min(tempdatosAmbos(:,2)); %Percentil 10
            end
            estadisPun(6,ipuntos)=sqrt(mean((tempdatosAmbos(:,2)-tempdatosAmbos(:,3)).^2)); %RMSE
            estadisPun(7,ipuntos)=estadisPun(6,ipuntos);                                            %NRMSE_A, it should be normalize with sigma2 in the excel template
            estadisPun(8,ipuntos)=sqrt(mean((tempdatosAmbos(:,2)-tempdatosAmbos(:,3)).^2./tempdatosAmbos(:,2).^2));   %NRMSE_B, normalized with each observation
            estadisPun(9,ipuntos)=estadisPun(6,ipuntos)/estadisPun(1,ipuntos);          %NRMSE1, normalized with observed ranges
            estadisPun(10,ipuntos)=estadisPun(6,ipuntos)/estadisPun(2,ipuntos);         %NRMSE2, normalized with observed means
            estadisPun(12,ipuntos)=mean(abs(tempdatosAmbos(:,2)-tempdatosAmbos(:,3)));  %MD, mean difference
            estadisPun(11,ipuntos)=estadisPun(6,ipuntos)/ estadisPun(12,ipuntos);                    %NRMSE_C, con mean difference
            estadisPun(13,ipuntos)=mean(abs((tempdatosAmbos(:,3)-tempdatosAmbos(:,2))./tempdatosAmbos(:,2).^2));      %NMD1=MAPE, normalized with each observation
            estadisPun(14,ipuntos)=estadisPun(12,ipuntos);                              %NMD2, it should be normalize with sigma2 in the excel template
            estadisPun(15,ipuntos)=estadisPun(12,ipuntos)/estadisPun(1,ipuntos);        %NMD3, normalized with observed ranges
            estadisPun(16,ipuntos)=estadisPun(12,ipuntos)/estadisPun(2,ipuntos);        %NMD4, normalized with observed means
            tempa=fitlm(tempdatosAmbos(:,2),tempdatosAmbos(:,3));
            estadisPun(17,ipuntos)=tempa.Rsquared.Ordinary;
            estadisPun(18,ipuntos)=tempa.Coefficients{2,1};
            estadisPun(19,ipuntos)=tempa.Coefficients{1,1};
            clear tempa
            
        end
    end
    %Saving the acum deformation series
    if ipuntos==1
        if exist('tempdatosAmbos','var')
            DefAcum=tempdatosAmbos;
        else
            DefAcum=zeros(1,3);
        end
    else
        if ~exist('tempdatosAmbos','var')
            tempdatosAmbos=zeros(1,3);
        end
        tempja=max(size(DefAcum,1),size(tempdatosAmbos,1));
        DefAcum=[DefAcum; zeros(tempja-size(DefAcum,1),3)];
        tempdatosAmbos2=[tempdatosAmbos; zeros(tempja-size(tempdatosAmbos,1),3)];
        DefAcum=[DefAcum tempdatosAmbos2];
    end
    clear tempdatosAmbos2;
    clear tempdatosAmbos
end

%% Velocity Statistics
estadisVel=nan(19,1);
if sum(sum(~isnan(veloParcial)))>0  %run only if there are points to analyze
    %assumes observation at source 1 and simulated at source 2
    P1= veloParcial (:,1);
    P2= veloParcial (:,2);
    
    [tempD,~]=find(isnan(P1));   %removing all nan
    P1(tempD,:)=[];
    [tempD,~]=find(isnan(P2));   %removing all nan
    P2(tempD,:)=[];
    
    MD= mean(abs(P2-P1));         %mean diference
    estadisVel(1,1)=max(P1)-min(P1);  %absolute range
    estadisVel(2,1)=abs(mean (P1));   %absolut mean
    estadisVel(3,1)=size(P1,1);       %nº observations
    if size(P1,1)>10
        estadisVel(4,1)= prctile(P1,90); %Percentil 90
        estadisVel(5,1)= prctile(P1,10); %Percentil 10
    else
        estadisVel(4,1)= max(P1); %Percentil 90
        estadisVel(5,1)= min(P1); %Percentil 10
    end
    estadisVel(6,1)=sqrt(mean((P1-P2).^2));   %RMSE
    estadisVel(7,1)=estadisVel(6); %NRMSE_A, it should be normalize with sigma2 in the excel template
    estadisVel(8,1)=sqrt(mean((P1-P2).^2./P1.^2));   %NRMSE_B, normalized with each observation
    estadisVel(9,1)=estadisVel(6)/estadisVel(1);  %NRMSE1, normalized with observed ranges
    estadisVel(10,1)=estadisVel(6)/estadisVel(2); %NRMSE2, normalized with observed means
    estadisVel(11,1)=estadisVel(6)/MD;            %NRMSE_C, normalized with mean difference
    estadisVel(12,1)=MD;                          %mean difference
    estadisVel(13,1)=mean(abs((P2-P1)./P1.^2));      %NMD1=MAPE, normalized with each observation
    estadisVel(14,1)=MD;                  %NMD2, it should be normalize with sigma2 in the excel template
    estadisVel(15,1)=MD/estadisVel(1);            %NMD3, normalized con with observed ranges
    estadisVel(16,1)=MD/estadisVel(2);            %NMD4, normalized with observed means
    tempa=fitlm(P1,P2);
    estadisVel(17,1)=tempa.Rsquared.Ordinary;
    estadisVel(18,1)=tempa.Coefficients{2,1};       %slope,
    estadisVel(19,1)=tempa.Coefficients{1,1};       %free term
    clear tempa
end

estadisAll=[estadisVel mean(estadisPun,2) estadisPun];


%% Saving
if sum(sum(~isnan(veloParcial)))>0
    %Creating the folder and file to collect the results
    carpetanueva=[fuente1 '_' fuente2];
    mkdir([pwd '\Results\' carpetanueva]);
    %             mkdir([pwd '\Results_' folder '\' carpetanueva]); -->-->-->para crear carpeta diferente en cada study area
    archivosalida=[pwd '\Results\' carpetanueva '\Validation_' fuente1 '_' fuente2 '.xlsx'];
    copyfile([pwd '\templateResults.xlsx'], archivosalida, 'f');
    
    %Saving velocity results
    tempescribir=[(1:npuntos)' veloParcial];
    writematrix(tempescribir, archivosalida,'Sheet','Velocity','Range','b3');
    tempescribir={'Loc' fuente1 fuente2};    %headings
    writecell(tempescribir, archivosalida,'Sheet','Velocity','Range','b2');
    clear tempescribir
    
    %Savind deformation time series
    tempescribir=[(1:npuntos)' veloParcial];
    writematrix(tempescribir, archivosalida,'Sheet','Velocity','Range','b3');
    for ipuntos=1:npuntos
        if ipuntos==1
            tempescribir={'Time' fuente1 fuente2};    %headings
            tempescribir2=ipuntos*ones(1,3);            %headings2
        else
            tempescribir=[tempescribir {'Time' fuente1 fuente2}];    %headings
            tempescribir2=[tempescribir2 ipuntos*ones(1,3)];        %headings2
        end
        
    end
    writecell(tempescribir, archivosalida,'Sheet','Deformation','Range','b3');
    writematrix(tempescribir2, archivosalida,'Sheet','Deformation','Range','b4');
    ja=DefAcum(:,1:3:3*npuntos);
    [tempf]=find(ja~=0);
    ja(tempf)=m2xdate(ja(tempf));
    DefAcum(:,1:3:3*npuntos)= ja;
    clear ja tempf
    writematrix(DefAcum, archivosalida,'Sheet','Deformation','Range','b5');
    clear tempescribir tempescribir2
    
    %Saving Statistics
    writematrix(estadisAll, archivosalida,'Sheet','Statis','Range','C5');
    
    
    %% Figures
    fig01=1; %Fig01 velocities
    if fig01==1
        figure (1)
        hgcf=gcf;
        hgcf.Position=[2200 380 600 400];
        hgca=gca;
        h1=plot(veloParcial(:,1),veloParcial(:,2));
        h1.LineStyle='none';
        h1.Marker='d';
        h1.MarkerFaceColor='k';
        h1.MarkerEdgeColor='k';
        
        hold on
        %1:1 line
        vmin=min(min(veloParcial));
        vmin2=vmin-1;
        h2=plot([1 vmin2],[1 vmin2]);
        h2.LineStyle='-';
        h2.Marker='none';
        h2.Color=[0.2 0.2 0.2];
        
        %Grey lines, five lines
        for i=1:5
            %                 h3(i)=plot([1 vmin2],[1 vmin2]+0.25*i);
            %                 h3(i).LineStyle=':';
            %                 h3(i).Marker='none';
            %                 h3(i).Color=[0.2*i 0.2*i 0.2*i];
            %
            %                 h4(i)=plot([1 vmin2],[1 vmin2]-0.25*i);
            %                 h4(i).LineStyle=':';
            %                 h4(i).Marker='none';
            %                 h4(i).Color=[0.2*i 0.2*i 0.2*i];
            %
        end
        
        %Grey lines, only one line
        h3=plot([1 vmin2],[1 vmin2]+1);
        h3.LineStyle=':';
        h3.Marker='none';
        h3.Color=[0.5 0.5 0.5];
        %Downlimit
        h4=plot([1 vmin2],[1 vmin2]-1);
        h4.LineStyle=':';
        h4.Marker='none';
        h4.Color=[0.5 0.5 0.5];
        
        % Adjustment line
        lineaajuste=1;
        if lineaajuste==1
            tempy=[vmin2 1]*estadisVel(18,1)+estadisVel(19,1);
            h5=plot([vmin2 1],tempy);
            h5.LineStyle='--';
            h5.Marker='none';
            h5.Color='r';
        end
        
        %Title of the axes
        hgca.XLabel.String=[fuente1 ' (cm/year)'];
        hgca.YLabel.String=[fuente2 ' (cm/year)'];
        
        %Limits of the axes
        hgca.XLim=[fix(vmin2) 1];
        hgca.YLim=[fix(vmin2) 1];
        
        linea=[pwd '\Results\' carpetanueva '\Fig01_Vel_' fuente1 '_' fuente2];
        print(gcf,linea,'-dtiff','-r600');
        close gcf
    end
    
    fig02=1;  %Fig02 Time series of each location
    if fig02==1
        for ipuntos=1:npuntos
            %     ipuntos=3;
            col=1+(ipuntos-1)*2;
            tempdatos1=datoscomp1(:,col:col+1);
            tempdatos2=datoscomp2(:,col:col+1);
            [tempf,~]=find(isnan(tempdatos1));   %removing all nan
            tempdatos1(tempf,:)=[];
            [tempf,~]=find(isnan(tempdatos2));   %removing all nan
            tempdatos2(tempf,:)=[];
            
            if ~isempty(tempdatos1) && ~isempty(tempdatos2)  %Avoid comparing points without data
                %Detection of contemporary data:
                if max(tempdatos1(1,1),tempdatos2(1,1))< min(tempdatos1(end,1),tempdatos2(end,1))
                    %Detection of the first begining:
                    if tempdatos1(1,1)<tempdatos2(1,1)
                        incremento2=interp1(tempdatos1(:,1),tempdatos1(:,2), tempdatos2(1,1));
                        tempdatos2modi=tempdatos2(:,2)+(incremento2-tempdatos2(1,2));
                        flag=1;
                    else
                        incremento1=interp1(tempdatos2(:,1),tempdatos2(:,2), tempdatos1(1,1));
                        tempdatos1modi=tempdatos1(:,2)+(incremento1-tempdatos1(1,2));
                        flag=2;
                    end
                end
                
                figure (1)
                hgcf=gcf;
                hgcf.Position=[2200 380 800 400];
                hgca=gca;
                %serie 1
                if flag==1
                    h1=plot(tempdatos1(:,1),tempdatos1(:,2));
                else
                    h1=plot(tempdatos1(:,1),tempdatos1modi(:,1));
                end
                h1.LineStyle='-';
                h1.Marker='d';
                h1.MarkerFaceColor='k';
                h1.MarkerEdgeColor='k';
                h1.MarkerSize=2;
                %serie 2
                hold on
                if flag==1
                    h2=plot(tempdatos2(:,1),tempdatos2modi(:,1));
                else
                    h2=plot(tempdatos2(:,1),tempdatos2(:,2));
                end
                h2.LineStyle='-';
                h2.Marker='d';
                h2.MarkerFaceColor='r';
                h2.MarkerEdgeColor='r';
                h2.MarkerSize=2;
                
                %Axes
                hgca.XLabel.String='Date';
                hgca.YLabel.String=['Deformation at location ' num2str(ipuntos) ' (cm)'];
                
                x3=[tempdatos1(:,1);tempdatos2(:,1)];
                x3=sort(x3);
                extension=0.03*(x3(end)-x3(1));
                hgca.XLim=[x3(1)-extension x3(end)+extension];
                hgca.XTick=datenum(year(x3(1)):year(x3(end)),1,1);
                hgca.XTickLabel=datestr(datenum(year(x3(1)):year(x3(end)),1,1),'yyyy');
                hgca.XLabel.String='Date';
                
                %Legend
                hl=legend([h1,h2],fuente1, fuente2);
                hl.FontSize=9;
                
                linea=[pwd '\Results\' carpetanueva '\Fig02_Def_' fuente1 '_' fuente2 '_Location ' num2str(ipuntos)];
                print(gcf,linea,'-dtiff','-r600');
                close gcf
                
            end
        end
    end
    
    fig03=1;  %insitu versus remote observations
    if fig03==1
        for ipuntos=1:npuntos
            datos=DefAcum(:,1+3*(ipuntos-1):3+3*(ipuntos-1));
            [tempf,~]=find(datos==0);   %removing all zeros
            datos(tempf,:)=[];
            [tempf,~]=find(isnan(datos));   %removing all zeros
            datos(tempf,:)=[];
            clear tempf
            if size(datos,1)>1  %at least two points to build the figure
                % figure for each location
                
                figure (1)
                hgcf=gcf;
                hgcf.Position=[2200 380 600 400];
                hgca=gca;
                h1=plot(datos(:,2),datos(:,3));
                h1.LineStyle='none';
                h1.Marker='d';
                h1.MarkerFaceColor='k';
                h1.MarkerEdgeColor='k';
                
                hold on
                %1:1 line
                xmin=min(min(datos(:,2:3)));
                xmax=max(max(datos(:,2:3)));
                xmin2=xmin-0.1*(xmax-xmin);
                xmax2=xmax+0.1*(xmax-xmin);
                h2=plot([xmin2 xmax2],[xmin2 xmax2]);
                h2.LineStyle='-';
                h2.Marker='none';
                h2.Color=[0.2 0.2 0.2];
                
                %Grey lines, only one line
                h3=plot([xmin2 xmax2],[xmin2 xmax2]+1);
                h3.LineStyle=':';
                h3.Marker='none';
                h3.Color=[0.5 0.5 0.5];
                %Downlimit
                h4=plot([xmin2 xmax2],[xmin2 xmax2]-1);
                h4.LineStyle=':';
                h4.Marker='none';
                h4.Color=[0.5 0.5 0.5];
                
                % adjustment line
                lineaajuste=1;
                if lineaajuste==1
                    tempy=[xmin2 xmax2]*estadisPun(18,ipuntos)+estadisPun(19,ipuntos);
                    h5=plot([xmin2 xmax2],tempy);
                    h5.LineStyle='--';
                    h5.Marker='none';
                    h5.Color='r';
                end
                
                %titles of the axes
                hgca.XLabel.String=[fuente1 ' (cm)'];
                hgca.YLabel.String=[fuente2 ' (cm)'];
                
                %Limits of the axes
                hgca.XLim=[xmin2 xmax2];
                hgca.YLim=[xmin2 xmax2];
                clear xmin xmax xmin2 xmax2
                %Saving figure
                linea=[pwd '\Results\' carpetanueva '\Fig03_CompareDef_' fuente1 '_' fuente2 '_Location ' num2str(ipuntos)];
                print(gcf,linea,'-dtiff','-r600');
                close gcf
            end
        end
    end
end