% loading data from the current period

% loading MNC and EOG data
disp('Loading MNC and EOG data')
for i=1:length(period)
    if exist(['mncdata/',stationname,datestr(period(i),'YY'),datestr(period(i),'mm'),datestr(period(i),'dd'),'.MNC'])>0
        MNCtemp=load(['mncdata/',stationname,datestr(period(i),'YY'),datestr(period(i),'mm'),datestr(period(i),'dd'),'.MNC']);
        MNCdata(i,:)=[period(i) MNCtemp(4) MNCtemp(6)];
    else
        MNCdata(i,:)=[period(i) -999.99 -999.99];
    end
    if exist(['eogdata/',stationname,datestr(period(i),'YY'),datestr(period(i),'mm'),datestr(period(i),'dd'),'.EOG'])>0
        EOGtemp=load(['eogdata/',stationname,datestr(period(i),'YY'),datestr(period(i),'mm'),datestr(period(i),'dd'),'.EOG']);
        EOGdata(i,:)=[period(i) EOGtemp];
    else
        EOGdata(i,:)=[period(i) -999.99];
        disp(['loading data: ',datestr(period(i),29)]);
    end    
end

% selecting which traj-data to load
for i=1:length(period)
    if EOGdata(i,2)>-999.99 & MNCdata(i,2)>-999.99
        disp(['loading trajectory data: ',datestr(period(i),29)]);
        % selecting which times we need to load back trajectories
        % starttime of formation of 1.7 nm diameter particles at Hyytiälä
        starttime(i)=MNCdata(i,2)-timeshift/24;
        % start time of formation of 1.7 nm particles at Hyytiälä average
        % to nearest full hour
        starttimefloor(i)=period(i)+floor((starttime(i)-floor(starttime(i)))*24)/24;
        % end time of formation of 1.7 nm diameter particles at Hyytiälä
        endtime(i)=MNCdata(i,3)-timeshift/24;
        % end time of formation of 1.7 nm particle at Hyytiälä average to
        % nearest full hour
        endtimeceil(i)=period(i)+ceil((endtime(i)-floor(starttime(i)))*24)/24;
        % time when the observed growth at Hyytiälä can no longer be
        % observed
        EOGtime(i)=EOGdata(i,2);
        % The first selected time when we start to calculate back
        % trajectories (in full hours, since trajectories can only be
        % calculated every full hour)
        starttraj(i)=period(i)+ceil((endtime(i)-floor(endtime(i))+1/24)*24)/24;
        % the last selected time for which we calculate back trajectories
        % (equal to EOGtime, but averaged to full hours)
        endtraj(i)=floor(EOGtime(i))+ceil((EOGtime(i)-floor(EOGtime(i)))*24)/24;
    else
        starttime(i)=-999.99;
        starttimefloor(i)=-999.99;
        endtime(i)=-999.99;
        endtimeceil(i)=-999.99;
        EOGtime(i)=-999.99;
        starttraj(i)=-999.99;
        endtraj(i)=-999.99;        
    end
end

% calculating trajdata for each day
latint3D=-999.99*ones(length(period),150,20);lonint3D=-999.99*ones(length(period),150,20);winddir3D=-999.99*ones(length(period),150,20);
for i=1:length(period)
    disp(['calculating trajectory data: ',datestr(period(i),29)]);
    % checking if we have event and if we have defined EOG time
    if starttraj(i)>-999.99
        trajdata=[];n=0;latint=[];lonint=[];latint2=[];lonint2=[];winddir=[];lattemp=[];lontemp=[];
        % from the start trajectory time to the end trajectory time
        for j=starttraj(i):1/24:endtraj(i)
            % index for each hour
            n=n+1;
            % loading the trajectory for each hour
            [slask1,slask2,slask3,slask4,day,hour,slask5,slask6,hourback,lat,long,z,p,ptemp,atemp,prec,mixz,RH]=...
                textread(['trajdata/',datestr(j+GMT_shift/24,'YYYY'),datestr(j+GMT_shift/24,'mm'),datestr(j+GMT_shift/24,'dd'),'_',datestr(j+GMT_shift/24,'HH'),'.dat'],...
                '%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%*[^\n]','headerlines',9);
            % how many hours backwards from the end of growth rate time should we go to reach the
            % start time of the 1.7 nm formation
            backtimemax=round((j-starttimefloor(i))*24);
            % how many hours backwards from the end of growth rate time
            % should we go to reach the end time of the 1.7 nm formation
            backtimemin=round((j-endtimeceil(i))*24);
            % this is the time interval (each hour) for the current
            % trajectory that defines where the 1.7 nm diaemter formation 
            % takes place
            timetemp=j-backtimemax/24:1/24:j-backtimemin/24;
            % this is the latitude intterval of where the formation takes
            % place for the current trajectory
                        
            %check if the trajectory data range is large enough (some
            %.dat files might be empty)
            if length(lat)>(backtimemax*3+1+3) %
                
                lattemp=lat(backtimemin*3+1:3:backtimemax*3+1);
                % this is the longitude intterval of where the formation takes
                % place for the current trajectory
                lontemp=long(backtimemin*3+1:3:backtimemax*3+1);
                % we interpolate the time interval in a finer mesh (each
                % 10 minutes)
                timeint=starttime(i):1/24:endtime(i);
                if abs(timeint(end)-endtime(i))>1/24/5
                    timeint(end+1)=endtime(i);
                end
                % then we get a new finer mesh for the latitute and longitude
                latint(n,:)=interp1(timetemp,lattemp,timeint);
                lonint(n,:)=interp1(timetemp,lontemp,timeint);
                % first we need to extend the trajectory with one starting and
                % one end point
                latint2(n,:)=interp1(j-backtimemax/24-1/24:1/24:j-backtimemin/24+1/24,lat(backtimemin*3+1-3:3:backtimemax*3+1+3),[timeint(1)-1/24 timeint timeint(end)+1/24]);
                lonint2(n,:)=interp1(j-backtimemax/24-1/24:1/24:j-backtimemin/24+1/24,long(backtimemin*3+1-3:3:backtimemax*3+1+3),[timeint(1)-1/24 timeint timeint(end)+1/24]);
                % then we calculate the wind direction for all the latint and
                % lonint points
                for k=1:size(latint,2);
                    [RNG,winddir(n,k)]=distance(latint2(n,k),lonint2(n,k),latint2(n,k+2),lonint2(n,k+2));
                end
            end
        end
       latint3D(i,1:size(latint,1),1:size(latint,2))=latint;
       lonint3D(i,1:size(latint,1),1:size(latint,2))=lonint;
       winddir3D(i,1:size(latint,1),1:size(latint,2))=winddir;
    end
end