% subroutine to calculate the grid count in each grid cell.

lower_lat_border=zeros(length(latlim(1):resolution_lat:latlim(2)-resolution_lat));
higher_lat_border=zeros(length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
lower_lon_border=zeros(length(lonlim(1):resolution_lon:lonlim(2)-resolution_lon));
higher_lon_border=zeros(length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));

gridvalues=zeros(length(lower_lat_border),length(lower_lon_border));
gridvalues_west=gridvalues;
gridvalues_north=gridvalues;
gridvalues_south=gridvalues;
gridvalues_east=gridvalues;

%temp grids
gridvaluestemp=gridvalues;
gridvalues_north_temp=gridvalues;
gridvalues_south_temp=gridvalues;
gridvalues_east_temp=gridvalues;
gridvalues_west_temp=gridvalues;

%% Calculate event-grid and event-winddirection
for k=1:length(period) % Loop over the full time period
        disp(['Calculating gridded version of new particle formation map: ', datestr(period(k),29)])
    for i=1:length(lower_lat_border);
        for j=1:length(lower_lon_border);
            if ~isempty(find((latint3D(k,:,:)>latlim(1)+i*resolution_lat & latint3D(k,:,:)<latlim(1)+i*resolution_lat+resolution_lat))); % & lonint3D>lonlim(1)+j*resolution_lon*j & lonint3D<lonlim(1)+resolution_lon*j+resolution_lon))
                indexi=find((latint3D(k,:,:)>latlim(1)+i*resolution_lat & latint3D(k,:,:)<latlim(1)+i*resolution_lat+resolution_lat));
                
                if ~isempty(find(lonint3D(k,indexi)>lonlim(1)+j*resolution_lon & lonint3D(k,indexi)<lonlim(1)+j*resolution_lon+resolution_lon));
                    indexij=indexi(find((lonint3D(k,indexi)>lonlim(1)+j*resolution_lon & lonint3D(k,indexi)<lonlim(1)+j*resolution_lon+resolution_lon)));
                    gridvaluestemp(i,j)=length(indexij);
                    
                    if ~isempty(winddir3D(k,indexij)>=225 & winddir3D(k,indexij)<315);
                        temp=(winddir3D(k,indexij)>=225 & winddir3D(k,indexij)<315);
                        gridvalues_west_temp(i,j)=length(find(temp>0));
                    end
                    if ~isempty(winddir3D(k,indexij)>=315 & winddir3D(k,indexij)<=360 | winddir3D(k,indexij)<45);
                        temp=(winddir3D(k,indexij)>=315 & winddir3D(k,indexij)<=360 | winddir3D(k,indexij)<45);
                        gridvalues_north_temp(i,j)=length(find(temp>0));
                    end
                    if ~isempty(winddir3D(k,indexij)>=135 & winddir3D(k,indexij)<225);
                        temp=(winddir3D(k,indexij)>=135 & winddir3D(k,indexij)<225);
                        gridvalues_south_temp(i,j)=length(find(temp>0));
                    end
                    if ~isempty(winddir3D(k,indexij)>=45 & winddir3D(k,indexij)<135);
                        temp=(winddir3D(k,indexij)>=45 & winddir3D(k,indexij)<135);
                        gridvalues_east_temp(i,j)=length(find(temp>0));
                    end
                end
            end
        end
    end
    
    % make all gridpoints with data>0 ==> =1 (this is done for each day
    % separately)
    for i2=1:size(gridvaluestemp,1)
        for i3=1:size(gridvaluestemp,2)
            if gridvaluestemp(i2,i3)>=1
                gridvaluestemp(i2,i3)=1;
            else
                gridvaluestemp(i2,i3)=0;
            end
            if gridvalues_east_temp(i2,i3)>=1
                gridvalues_east_temp(i2,i3)=1;
            else
                gridvalues_east_temp(i2,i3)=0;
            end
            if gridvalues_west_temp(i2,i3)>=1
                gridvalues_west_temp(i2,i3)=1;
            else
                gridvalues_west_temp(i2,i3)=0;
            end
            if gridvalues_south_temp(i2,i3)>=1
                gridvalues_south_temp(i2,i3)=1;
            else
                gridvalues_south_temp(i2,i3)=0;
            end
            if gridvalues_north_temp(i2,i3)>=1
                gridvalues_north_temp(i2,i3)=1;
            else
                gridvalues_north_temp(i2,i3)=0;
            end
        end
    end
    % Add gridvalues (now equal to 0 or 1)
    gridvalues=gridvalues+gridvaluestemp;
    gridvalues_north=gridvalues_north+gridvalues_north_temp;
    gridvalues_south=gridvalues_south+gridvalues_south_temp;
    gridvalues_east=gridvalues_east+gridvalues_east_temp;
    gridvalues_west=gridvalues_west+gridvalues_west_temp;
    gridvaluestemp=zeros(size(gridvaluestemp));
    gridvalues_west_temp=zeros(size(gridvalues_west_temp));
    gridvalues_north_temp=zeros(size(gridvalues_north_temp));
    gridvalues_east_temp=zeros(size(gridvalues_east_temp));
    gridvalues_south_temp=zeros(size(gridvalues_south_temp));
end

%%%%%%%

% figure of grid count in all wind directions
disp('Plotting figure 2')
figure(2);clf;
hold on
worldmap([latlim(1) latlim(2)],[lonlim(1) lonlim(2)]);
latmatrix=(latlim(1)+resolution_lat/2:resolution_lat:latlim(2)-resolution_lat/2)'*ones(1,length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));
lonmatrix=(lonlim(1)+resolution_lon/2:resolution_lon:lonlim(2)-resolution_lat/2)'*ones(1,length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0 0 1]);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.9 0.9 0.9]);
grid_log=gridvalues;
for i=1:size(latmatrix,1)
    for j=1:size(latmatrix,2);
        if grid_log(i,j)<=0
            grid_log(i,j)=-1;
        else
            grid_log(i,j)=log10(gridvalues(i,j));
        end
    end            
end
surfm(latmatrix,lonmatrix',grid_log);
% color bar
caxis([-0.7 max(max(grid_log))]);
m=[];m=colormap('jet');
m(1,:)=[0 0 0.25];
n=0;l=[];
for i=1:5:size(m,1);
    n=n+1;
    l(n,:)=m(i,:);
end
colormap(l);
l=colorbar('vert');
set(l,'position',[0.90 0.25 .03 .63])
set(l,'ytick',[0 1 2 3]);
set(l,'yticklabel',[1 10 100 1000]);
% denoting station with cross
linem(poslat,poslon,'X','linewidth',2,'markersize',20,'color',[0.99 0.99 0.99]);
%Title - station name and period
timestring=[datestr(min(period)) {'to'} datestr(max(period))];
titlestring=[num2str(stationname) {} timestring];
title(titlestring,'FontSize',12,'fontweight','bold')
distkm=500; %Distance of radius in km
distdeg = distdim(distkm,'kilometers','deg');
[latc,lonc] = scircle1(poslat,poslon,distdeg);
geoshow(latc,lonc,'color',[0.99 0.99 0.99]);
linem(poslat,poslon,'wX','linewidth',2,'markersize',20);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.99 0.99 0.99]);
hold off
print -djpeg -r300 map2
print -dmeta -r300 map2
saveas(figure(2),'map2.fig');

% figure of grid count of westerly winds
disp('Plotting figure 3')
figure(3);clf;
hold on
worldmap([latlim(1) latlim(2)],[lonlim(1) lonlim(2)]);
latmatrix=(latlim(1)+resolution_lat/2:resolution_lat:latlim(2)-resolution_lat/2)'*ones(1,length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));
lonmatrix=(lonlim(1)+resolution_lon/2:resolution_lon:lonlim(2)-resolution_lat/2)'*ones(1,length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0 0 1]);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.99 0.99 0.99]);
grid_log_west=gridvalues_west;
for i=1:size(latmatrix,1)
    for j=1:size(latmatrix,2);
        if grid_log_west(i,j)<=0
            grid_log_west(i,j)=-1;
        else
            grid_log_west(i,j)=log10(gridvalues_west(i,j));
        end
    end            
end
surfm(latmatrix,lonmatrix',grid_log_west);
% color bar
l=colorbar('vert');
set(l,'position',[0.90 0.25 .03 .63])
set(l,'ytick',[-1 0 1 2 3]);
set(l,'yticklabel',[0 1 10 100 1000]);
% denoting station with cross
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
%Title - station name and period
timestring=[datestr(min(period)) {'to'} datestr(max(period))];
titlestring=[num2str(stationname) {} timestring];
title(titlestring,'FontSize',12,'fontweight','bold')
distkm=500; %Distance of radius in km
distdeg = distdim(distkm,'kilometers','deg');
[latc,lonc] = scircle1(poslat,poslon,distdeg);
geoshow(latc,lonc,'color',[0.99 0.99 0.99]);
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
hold off
print -djpeg -r300 map3

% figure of grid count of northerly winds
disp('Plotting figure 4')
figure(4);clf;
hold on
worldmap([latlim(1) latlim(2)],[lonlim(1) lonlim(2)]);
latmatrix=(latlim(1)+resolution_lat/2:resolution_lat:latlim(2)-resolution_lat/2)'*ones(1,length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));
lonmatrix=(lonlim(1)+resolution_lon/2:resolution_lon:lonlim(2)-resolution_lat/2)'*ones(1,length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0 0 1]);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.99 0.99 0.99]);
grid_log_north=gridvalues_north;
for i=1:size(latmatrix,1)
    for j=1:size(latmatrix,2);
        if grid_log_north(i,j)<=0
            grid_log_north(i,j)=-1;
        else
            grid_log_north(i,j)=log10(gridvalues_north(i,j));
        end
    end            
end
surfm(latmatrix,lonmatrix',grid_log_north);
% color bar
l=colorbar('vert');
set(l,'position',[0.90 0.25 .03 .63])
set(l,'ytick',[-1 0 1 2 3]);
set(l,'yticklabel',[0 1 10 100 1000]);
% denoting station with cross
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
%Title - station name and period
timestring=[datestr(min(period)) {'to'} datestr(max(period))];
titlestring=[num2str(stationname) {} timestring];
title(titlestring,'FontSize',12,'fontweight','bold')
distkm=500; %Distance of radius in km
distdeg = distdim(distkm,'kilometers','deg');
[latc,lonc] = scircle1(poslat,poslon,distdeg);
geoshow(latc,lonc,'color',[0.99 0.99 0.99]);
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
hold off
print -djpeg -r300 map4

% figure of grid count of easterly winds
disp('Plotting figure 5')
figure(5);clf;
hold on
worldmap([latlim(1) latlim(2)],[lonlim(1) lonlim(2)]);
latmatrix=(latlim(1)+resolution_lat/2:resolution_lat:latlim(2)-resolution_lat/2)'*ones(1,length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));
lonmatrix=(lonlim(1)+resolution_lon/2:resolution_lon:lonlim(2)-resolution_lat/2)'*ones(1,length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0 0 1]);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.99 0.99 0.99]);
grid_log_east=gridvalues_east;
for i=1:size(latmatrix,1)
    for j=1:size(latmatrix,2);
        if grid_log_east(i,j)<=0
            grid_log_east(i,j)=-1;
        else
            grid_log_east(i,j)=log10(gridvalues_east(i,j));
        end
    end            
end
surfm(latmatrix,lonmatrix',grid_log_east);
% color bar
l=colorbar('vert');
set(l,'position',[0.90 0.25 .03 .63])
set(l,'ytick',[-1 0 1 2 3]);
set(l,'yticklabel',[0 1 10 100 1000]);
% denoting station with cross
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
%Title - station name and period
timestring=[datestr(min(period)) {'to'} datestr(max(period))];
titlestring=[num2str(stationname) {} timestring];
title(titlestring,'FontSize',12,'fontweight','bold')
distkm=500; %Distance of radius in km
distdeg = distdim(distkm,'kilometers','deg');
[latc,lonc] = scircle1(poslat,poslon,distdeg);
geoshow(latc,lonc,'color',[0.99 0.99 0.99]);
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
hold off
print -djpeg -r300 map5

% figure of grid count of southerly winds
disp('Plotting figure 6')
figure(6);clf;
hold on
worldmap([latlim(1) latlim(2)],[lonlim(1) lonlim(2)]);
latmatrix=(latlim(1)+resolution_lat/2:resolution_lat:latlim(2)-resolution_lat/2)'*ones(1,length(lonlim(1)+resolution_lon:resolution_lon:lonlim(2)));
lonmatrix=(lonlim(1)+resolution_lon/2:resolution_lon:lonlim(2)-resolution_lat/2)'*ones(1,length(latlim(1)+resolution_lat:resolution_lat:latlim(2)));
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0 0 1]);
geoshow([land.Lat], [land.Lon], 'Color' ,[0.99 0.99 0.99]);
grid_log_south=gridvalues_south;
for i=1:size(latmatrix,1)
    for j=1:size(latmatrix,2);
        if grid_log_south(i,j)<=0
            grid_log_south(i,j)=-1;
        else
            grid_log_south(i,j)=log10(gridvalues_south(i,j));
        end
    end            
end
surfm(latmatrix,lonmatrix',grid_log_south);
% color bar
l=colorbar('vert');
set(l,'position',[0.90 0.25 .03 .63])
set(l,'ytick',[-1 0 1 2 3]);
set(l,'yticklabel',[0 1 10 100 1000]);
% denoting station with cross
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
%Title - station name and period
timestring=[datestr(min(period)) {'to'} datestr(max(period))];
titlestring=[num2str(stationname) {} timestring];
title(titlestring,'FontSize',12,'fontweight','bold')
distkm=500; %Distance of radius in km
distdeg = distdim(distkm,'kilometers','deg');
[latc,lonc] = scircle1(poslat,poslon,distdeg);
geoshow(latc,lonc,'color',[0.99 0.99 0.99]);
linem(poslat,poslon,'kX','linewidth',2,'markersize',20);
hold off
print -djpeg -r300 map6

% saving data
disp('Saving data')
save('gridvalues.GRI','gridvalues','-double','-ascii');
save('gridvalues_west.GRI','gridvalues_west','-double','-ascii');
save('gridvalues_north.GRI','gridvalues_north','-double','-ascii');
save('gridvalues_south.GRI','gridvalues_south','-double','-ascii');
save('gridvalues_east.GRI','gridvalues_east','-double','-ascii');