clear;
clc;
%Sorter
%Current Configuration: Mexico Project


%Program Test Successfully
%If plot is x=date, y=altitude, z=altitude
%DATAREF = column 2, DATAGRAB = column 2
%Data will show 1:1 slope in y-z plane
%This shows that data does map ONLY available
%altitude with its corresponding data measurement




%Files Reference No.
%Ex: P0102021.WIN = P010 + Ref No. + .WIN
%                           Ref No. = 2021

%Mode: 0 = split data, 1 = read data into Matlab
Mode=0;
DATAREF=2;     %Reference data for y-axis
DATAGRAB=7;    %Plot data from this column
Ref_No=[2021;
    2032;
    2041;
    2051;
    2061;
    2062;
    2072;
    2082;
    2091;
    2092;
    2101;
    2102;
    2112;
    2121;
    2122;
    2131;
    2142;
    2151;
    2161;
    2162;
    2171;
    2172;
    2181;
    2182;
    2221;
    2222;
    2231;
    2232;
    2241;
    2242;
    2251;
    2261;
    2262;
    2271;
    2272;
    2281;
    2282;
    2291];

%Start of Code----------------------------------

%Find number of files
sizefactor=0;
max_filecount=size(Ref_No);

%Populate directory w/ data in correct format and load into Matlab
for filecount=1:1:max_filecount(1,1)
    filepath=strcat('P010',num2str(Ref_No(filecount)),'.WIN');

    %Splitting data in file
    pfile1=strcat('P010',num2str(Ref_No(filecount)),'_h.WIN');
    pfile2=strcat('P010',num2str(Ref_No(filecount)),'_d.WIN');

    if (Mode == 0)
        fid=fopen(filepath, 'r');
        fheaderid = fopen(pfile1,'w');
        fdataid = fopen(pfile2,'w');
        
        tline = fgets(fid);
        fprintf(fheaderid,tline);
        tline = fgets(fid);
        fprintf(fheaderid,tline);
        fclose(fheaderid);

        while feof(fid) == 0
            tline = fgets(fid);
            fprintf(fdataid,tline);
        end
        fclose(fdataid);
        fclose(fid);
    end
    %Load data
    output(filecount).data=load(pfile2);

    currentdata=output(filecount).data;
    currentsize=size(currentdata);

    if (currentsize(1,1)>sizefactor)
        sizefactor=currentsize(1,1);
    end

end

%Find max intervals
    datemax=max_filecount(1,1);
    arraymax=sizefactor;

%Populate meshgrid for contour plot
    [X, Y]= meshgrid(1:1:datemax,1:1:arraymax);
    Z=X;
    Z(:,:)=0;

%Merge values into mesh

for date_count=1:datemax
    information=output(date_count).data;
    
    sizeDATREF=size(information(:,DATAREF));
    inputA=information(:,DATAREF);
    if sizeDATREF(1,1)<sizefactor
        inputA((sizeDATREF(1,1)+1):sizefactor,:)=NaN;
    end
    
    sizeDATGRA=size(information(:,DATAGRAB));
    inputB=information(:,DATAGRAB);
    if sizeDATGRA(1,1)<sizefactor
        inputB((sizeDATGRA(1,1)+1):sizefactor,:)=NaN;
    end
    
        
    Y(:,date_count)=inputA;
    Z(:,date_count)=inputB;

end

for array_count=1:arraymax
    X(array_count,:)=Ref_No';
end

figure(1)
contour(X,Y,Z)
figure(2)
surf(X,Y,Z)