[points]=hdfreader_quikscat()

function [points]=hdfreader_quikscat()
% 
% SYNTAX: STRUCT_MATRIX=hdfreader_quikscat();
%
% This code reads as many files as you want (1 or >1) of your 
% choice from dataset PODAAC JPL QuikSCAT Seawinds L3 passings, 
% generating a struct-organized maxtrix.
%
% $$ FILES IN '*.gz' MUST BE EXTRACTED.
%
% $$ IS NOT NECESSARY TO RENAME THE EXTRACTED FILES.
%
% $$ On days when some point of the grid sample was not recovered its value,
% is placed NaN. (believe me, this will greatly facilitate the analysis).
%
% This code is based on PC Windows Platform, tested in Windows Seven
% Ultimate Edition (64-bit) on Matlab 64-bit version 7.8.0 (2009a) and
% same version of Windows 32-bit.
%
%
% 16/11/2009 Gustavo Hattenhauer Gomes ([e-mail])
% Student in Oceanography. Universidade do Vale do Itaja - Itaja/SC -
% Brazil.
%
% 20/11/2009 - add capability of analyse data of different years and
% without are in sequence.
%  ------------------------------------------------------------------

clear all;clc
% Specify GEOGRAPHIC COORDINATES for study area
% Longitude [0.125 to 359.875] ... Latitude [-89.875 to 89.875]
lat1= -25.625; lat2= -29.625; long1=313.125; long2=310.875;

% Create grid-box and maximun/mininum rows & colums
lat=(-89.875:0.25:89.875)'; long=(0.125:0.25:359.875)';
fa1=find(lat==lat1); fa2=find(lat==lat2); ta1=fa1; ta2=fa2;
fo1=find(long==long1); fo2=find(long==long2); to1=fo1; to2=fo2;

% Make sure that lat1 is greater lat2
if fa2 > fa1
    fa1=ta2;
    fa2=ta1;
end
% Make sure that long1 is greater long2
if fo1 > fo2
    fo1=to2;
    fo2=to1;
end

% Specify files to analyse
[filename,pathname]=uigetfile('*.*','MultiSelect','on');
if ischar(filename)==1
    filename.filename=filename;
    filename=struct2cell(filename);
end
for g=1:length(filename)
    yt=filename{g};
    y0=str2num(yt(11:14));
    yy(g,1)=y0;
end
y_min=min(yy); y_max=max(yy); year=y_min; f_length=length(filename);

% 'Last day'
endday=filename{length(filename)};

% start read algorithm with hdf_sd
X=0;
while X<=1
    year=y_min;
    if X==0
        PS='asc'; XX=0; XY=14;        
    elseif X==1
        PS='des'; XX=1; XY=15;
    end
    p1=[PS,'_avg_wind_speed']; p2=[PS,'_avg_wind_vel_u'];
    p3=[PS,'_avg_wind_vel_v']; p4=[PS,'_avg_wind_speed_sq'];
    p5=[PS,'_wvc_count']; p6=[PS,'_time_frac'];
    p7=[PS,'_rain_prob']; p8=[PS,'_rain_flag'];
    parameters.P1=p1; parameters.P2=p2; parameters.P3=p3;
    parameters.P4=p4; parameters.P5=p5; parameters.P6=p6;
    parameters.P7=p7; parameters.P8=p8;
    
    % start analyse
    day=1;
    while year <= y_max
        P=1;
        while day <= f_length
            FF=filename{day}; FD=str2num(FF(15:17)); FY=str2num(FF(11:14));
            if FY ~= year
                year=FY;
            end
            nday=filename{day};            
            pmts.data1=hdfread(nday,p1); pmts.data2=hdfread(nday,p2);
            pmts.data3=hdfread(nday,p3); pmts.data4=hdfread(nday,p4);
            pmts.data5=hdfread(nday,p5); pmts.data6=hdfread(nday,p6);
            pmts.data7=hdfread(nday,p7); pmts.data8=hdfread(nday,p8);
            
            u=1;
            for parameter=XX:2:XY
                dt5=pmts.data5; % asc or des wvc_count (is true or no)
                sd_id=hdfsd('start',nday,'read'); %open file HDF
                attr_idx1=hdfsd('findattr',sd_id,'date_of_average');
                [date_ave,status]=hdfsd('readattr',sd_id,attr_idx1);
                sds_idx=parameter; % relation with XX and XY
                sds_id=hdfsd('select',sd_id,sds_idx);
                [dsname, dsndims, dsdims,dstype,dsatts,stat]=hdfsd('getinfo',sds_id);
                attr_idx=parameter; % relation with XX and XY
                [scale_factor,status]=hdfsd('readattr',sds_id,0);% scale factor
                [add_offset,status]=hdfsd('readattr',sds_id,2);% offset
                status=hdfsd('endaccess',sds_id); % close acess to dataset
                status=hdfsd('end',sd_id); % close acess to hdf file
                
                dat=['data',num2str(u)];
                pmt=['P',num2str(u)];
                pmt=parameters.(pmt);
                data=pmts.(dat);
                data=double(data)*double(scale_factor) - double(add_offset);
                
                % Creating matrix with data         
                z=1;
                for W=fo1:fo2
                    for Y=fa2:fa1
                        if dt5(Y,W) > 0
                            new(z,1)=long(W,1);
                            new(z,2)=lat(Y,1);
                            new(z,u+2)=data(Y,W);
                            z=z+1;
                        elseif dt5(Y,W)==0
                            new(z,1)=long(W,1);
                            new(z,2)=lat(Y,1);
                            new(z,u+2)=NaN;
                            z=z+1;
                        end
                    end
                end                
                clear data
                u=u+1;
            end
            
            % Saving in structured matrix
            dd=['day',num2str(FD)]; aa=['year',num2str(year)];
            jj=[num2str(FD),'_',num2str(year),'_',PS];
            points.(aa).(dd).(PS)=new;
            %
            clear pmts new
            disp(jj); P=P+1;      
            if f_length == day
                year=year+100;
            end
            day=day+1;
        end
    end
    
    % Define what's the passing to analyse. Interspersing between asc or des
    if X==0
        X=1;
    elseif X==1
        X=0;
        if strcmp(nday,endday)==1
            X=2;
        end        
    end
    clear parameters
end
end