%%calculated for heading based on mean ground direction at departure 202.8
clear 
close all
%year=2014;
year=2017;
IsRadar=0; %0= mean heading 208.4 (red),n-1:GD=202.8; (green) n:HD=208.4 (blue)
           %1= variable ground direction taken from daily values radar; n-1=HD=208.4 (blue) n:GD=202.8 (green) 
           %2=mean ground direction at departure 202.8; n=mean head 208.4 (blue); n-1:GD=202.8 (green)  

mhead=208.4; % HD, mean heading from radar
mgroundD=202.8; %GD, mean ground direction from radar
           
           
%check on which daysradar was available % this part is left our as it was finally not used

load(['MeteoMatrix',num2str(year), '_925']); % load meteorological data (ERA40 at 925mbar pressure level)
load Coast % load shapefile coast
load Coast2 % load 2nd shapefile coast
nDays=62; % number of days of migration
nTracks=10000; % number of bird tracks
nSteps=60; % number of calculation steps (time integrations)
dt=0.5;%timestep in hours
if year == 2008 || year == 2012 || year == 2016
    DayS=228; %startday 15 Aug
else
    DayS=227; %startday 15 Aug
end
as=10;
Results=[];
TDist=0;
%LongDistTable at Allegranza
TabLong=[-20 -15:0.5:-10 -5];
TabDist=[-628.6 -144.3 -95.9 -47.5 1 49.4 97.8 146.3 194.7 243.2 291.6 340 824.2];
longlat=[];
% % image width and height
imgWidthPix = 1900;
imgHeightPix = 1050;
% set the image sizes, resolution and units; also see print option -r100
imgResolution = 100; % 100 dpi
imgWidth = imgWidthPix / imgResolution;
imgHeight = imgHeightPix / imgResolution;
% this line is necessary for rendering without openGL drivers/physical screen
set(0, 'DefaultFigureRenderer', 'zbuffer');

% for suncycle interpolation
SSLong=(-9.7:1:-4.7); SSLong(2,:)=SSLong(1,:);
SSLat=[36.5 36.5 36.5 36.5 36.5 36.5 ; 38 38 38 38 38 38];

%initiation matrices
Time(1:nTracks,1:nSteps,1:nDays)=NaN;
Lat(1:nTracks,1:nSteps,1:nDays)=NaN;
Long(1:nTracks,1:nSteps,1:nDays)=NaN;
wu(1:nTracks,1:nSteps,1:nDays)=NaN;
wv(1:nTracks,1:nSteps,1:nDays)=NaN;
gu(1:nTracks,1:nSteps,1:nDays)=NaN;
gv(1:nTracks,1:nSteps,1:nDays)=NaN;
au(1:nTracks,1:nSteps,1:nDays)=NaN;
av(1:nTracks,1:nSteps,1:nDays)=NaN;
GS(1:nTracks,1:nSteps,1:nDays)=NaN;
TDist(1:nTracks,1:nSteps,1:nDays)=NaN;
HD(1:nTracks,1:nSteps,1:nDays)=NaN;
Distu_unit(1:nTracks,1:nSteps,1:nDays)=NaN;
Distv_unit(1:nTracks,1:nSteps,1:nDays)=NaN;



scrsz = get(0,'ScreenSize');
Daynr=43; % for figure 2 only
for ii=Daynr:Daynr %only run for one day otherwise 1:nDays for nDays
    Tdum=DayS-2+ii+datenum([num2str(year),'-01-01 00:00:00']);
    TT=datestr(Tdum, 'yyyy mm dd'); TT(5)=[]; TT(7)=[];
 
    Long(1:nTracks,1,ii)=-8.7+3*rand(1,nTracks)';Long(nTracks-1:nTracks, 1, ii)=-7.9517;
    Lat(1:nTracks,1,ii)=37 +0.5*rand(1,nTracks)'; Lat(nTracks-1:nTracks, 1, ii)=37.3053;
    SSTime(2,6)=0;
    for i=1:nSteps
        if i==1
            for ic=1:6
                for ir=1:2  % suncycle starting points
                   DUM(1:2)=suncycle(SSLat(ir,ic), SSLong(ir,ic), [str2double(TT(1:4)), str2double(TT(5:6)), str2double(TT(7:8))], 2880)/24; 
                   SSTime(ir,ic)=DUM(2);
                end
            end 
            for k=1:nTracks
                  %calculate for each track the sunset time based on LatLong
                 DumTime=interpn(SSLat,SSLong,SSTime,Lat(k,1,ii),Long(k,1,ii),'linear', -999);
                 Time(k,1,ii)=DayS-1+ii+DumTime;
            end
            meanSTime(ii)= mean(Time(:,1,ii));
        else
           Time(1:nTracks,i,ii)=Time(1:nTracks,i-1,ii)+0.5/24;
        end
        
        %interpolate wu and wv 
        wu(1:nTracks,i,ii)=interpn(Mlat/100,Mlong/100,MTime, Mu, Lat(1:nTracks,i,ii),Long(1:nTracks,i,ii),Time(1:nTracks,i,ii),'linear', -999);
        wv(1:nTracks,i,ii)=interpn(Mlat/100,Mlong/100,MTime, Mv, Lat(1:nTracks,i,ii),Long(1:nTracks,i,ii),Time(1:nTracks,i,ii),'linear', -999);
        
        if i==1
           HD(1:nTracks-2,1,ii)=mhead+5*randn(nTracks-2,1);
           f2=(90:5:270); % heading (orientation in air)
           f1=(as*sin(f2/360*2*pi)+wu(nTracks-1,i,ii))./(as*cos(f2/360*2*pi)+wv(nTracks-1,i,ii)); %tan(GD)
           HD(nTracks-1,1,ii)=interp1(f1,f2,0.42036);
           HD(nTracks,1,ii)=mhead;
        end

        au(1:nTracks,i,ii)=as*sin(HD(1:nTracks,1,ii)/360*2*pi);
        av(1:nTracks,i,ii)=as*cos(HD(1:nTracks,1,ii)/360*2*pi); 
        gu(1:nTracks,i,ii)=au(1:nTracks,i,ii)+wu(1:nTracks,i,ii); % groundspeed u
        gv(1:nTracks,i,ii)=av(1:nTracks,i,ii)+wv(1:nTracks,i,ii); % groundspeed v
        GS(1:nTracks,i,ii)=(gu(1:nTracks,i,ii).*gu(1:nTracks,i,ii)+gv(1:nTracks,i,ii).*gv(1:nTracks,i,ii)).^0.5; %ground speed [m/s]
      
        if i==1
           TDist(1:nTracks,i,ii)=GS(1:nTracks,i,ii)*1800/1000;  
        else
           TDist(1:nTracks,i,ii)=TDist(1:nTracks,i-1,ii)+GS(1:nTracks,i,ii)*1800/1000; % total route distance [km]
        end

        for j=1:nTracks
           Distu_unit(j,i,ii)=distWB([Lat(j,i,ii) Long(j,i,ii)-0.5] ,[Lat(j,i,ii) Long(j,i,ii)+0.5]); %km per degree
           Distv_unit(j,i,ii)=distWB([Lat(j,i,ii)-0.5 Long(j,i,ii)] ,[Lat(j,i,ii)+0.5 Long(j,i,ii)]); %km per degree
        end
        if i<nSteps
            Lat(1:nTracks,i+1,ii)=Lat(1:nTracks,i,ii)+gv(1:nTracks,i,ii).*1800./1000./Distv_unit(1:nTracks,i,ii);
            Long(1:nTracks,i+1,ii)=Long(1:nTracks,i,ii)+gu(1:nTracks,i,ii).*1800./1000./Distu_unit(1:nTracks,i,ii);
        end
    end%for i=1:60 time steps
 
    WU=wu(:,:,ii); LAT=Lat(:,:,ii); LONG=Long(:,:,ii); TIME=Time(:,:,ii);
    dum=find(WU<-900); LAT(dum)=NaN;LONG(dum)=NaN; TIME(dum)=NaN;
    dum2=(LAT>29.4); dum3=sum(dum2')';
    wu(:,:,ii)=WU; Lat(:,:,ii)=LAT; Long(:,:,ii)=LONG; Time(:,:,ii)=TIME;
   
    
    for n=1:nTracks
       if dum3(n)<60 && dum3(n)>0
           ArrLong(n,1,ii)=((29.4-Lat(n,dum3(n),ii))./(Lat(n,dum3(n)+1,ii)-Lat(n,dum3(n),ii))).*(Long(n,dum3(n)+1,ii)-Long(n,dum3(n),ii))+ Long(n,dum3(n),ii);
           ArrTime(n,1,ii)=((29.4-Lat(n,dum3(n),ii))./(Lat(n,dum3(n)+1,ii)-Lat(n,dum3(n),ii))).*(Time(n,dum3(n)+1,ii)-Time(n,dum3(n),ii))+ Time(n,dum3(n),ii);
       else
           ArrLong(n,1,ii)=NaN;
           ArrTime(n,1,ii)=NaN;
       end
    end
    ArrDist(:,1,ii)=interp1(TabLong, TabDist,ArrLong(:,1,ii));
    WB1=ArrLong(1:nTracks-2,1,ii)<-13.0 ; WB2=ArrLong(1:nTracks-2,1,ii)>-14.0; Allegranza=find(WB1.*WB2 ==1);% arrived on allegranza
    if isempty(Allegranza)
        birddens(ii)=0; 
        AT(ii)=NaN;
        MedianDist(ii)=nanmedian(ArrDist(1:nTracks-2,1,ii));
        MedianDistA(ii)=NaN;
    else
        birddens(ii)=floor(length(Allegranza)/nTracks*1000)/10;
        ALong=ArrLong(Allegranza,1,ii);
        ATime=ArrTime(Allegranza,1,ii);
        AT(ii)=(mean(ATime)-meanSTime(ii))*24;
        MedianDistA(ii)=median(ArrDist(Allegranza,1,ii));
        MedianDist(ii)=nanmedian(ArrDist(1:nTracks-2,1,ii));
    end
    Dist5(ii)=prctile(ArrDist(1:nTracks-2,1,ii),5)';
    Dist95(ii)=prctile(ArrDist(1:nTracks-2,1,ii),95)';    
 
%% the making of plots for 62 days is here commented out    
%     %make subplot
%     if ii==1 || ii==33
%       scrsz = get(0,'ScreenSize');
%       figure('Position',[30 scrsz(4)/20 scrsz(3)-50 scrsz(4)/1.2]);      
%     end
%     if ii>32
%        subplot(4,8,ii-32)
%     else
%        subplot(4,8,ii)
%     end
    
    %draw map
%     for i=1:7
%        coastlon=Coast(i).X; coastlat=Coast(i).Y;
%        plot(coastlon, coastlat, '.k-'); 
%        hold on
%     end
%     for i=1:15
%        coastlon=Coast2(i).X; coastlat=Coast2(i).Y;
%        plot(coastlon, coastlat, '.k-'); 
%        hold on
%     end
%     xlim([-15 -5]); ylim([28 38])
%         for j=1:40
%             h3=plot(Long(j,1:60,ii),Lat(j,1:60,ii),'-','color',[0.5 0.5 0.5]);
%             h3.Color(4) = 0.1;
%             hold on
%         end
% 
%         if ii==1
%             title([TT,', P:',F2(17:19),' R:',num2str(IsRadar)],'fontsize',9)
%         else
%             if IsRadar==0
%                title([TT, ', AS',num2str(as),' HD=208'],'fontsize',9)
%             elseif IsRadar==1
%                title([TT, ', MT',num2str(MT(ii)),' GD=',num2str(round(GD(nTracks-1,1,ii))) ],'fontsize',9)
%             elseif IsRadar==2
%                title([TT, ', AS',num2str(as),' GD=203'],'fontsize',9)
% 
%             end
%         end
%         xlabel(['B',num2str(birddens(ii)),';T',num2str(round(AT(ii)*100)/100),';M',num2str(round(MedianDist(ii)*10)/10), ...
%                 ';',num2str(round(MedianDistA(ii)*10)/10),';',num2str(round(Dist5(ii)*10)/10)],'fontsize',9);
        Results=[Results; ii-1+DayS birddens(ii) AT(ii) MedianDist(ii) MedianDistA(ii) Dist5(ii) Dist95(ii)];

%         drawnow
   
        
%         if ii==32
%            if IsRadar==0
%               Filename=['D0_',TT(1:4),'-',F2(12:15),'-1-',num2str(as)];
%            elseif IsRadar==1
%               Filename=['D1_',TT(1:4),'-',F2(12:15),'-1-',num2str(as)];
%            elseif IsRadar==2
%                Filename=['D2_',TT(1:4),'-',F2(12:15),'-1-',num2str(as)];
%            end
%              saveas(gcf,[Filename, '.fig'])
%             set(gcf, 'PaperUnits', 'inches');
%             set(gcf, 'PaperPosition', [0 0 imgWidth imgHeight]);
%             print ('-dpng', '-r100', Filename)
%          elseif ii==62 
%             if IsRadar==0
%               Filename=['D0_',TT(1:4),'-',F2(12:15),'-2-',num2str(as)];
%             elseif IsRadar==1
%               Filename=['D1_',TT(1:4),'-',F2(12:15),'-2-',num2str(as)];
%             elseif IsRadar==2
%                Filename=['D2_',TT(1:4),'-',F2(12:15),'-2-',num2str(as)];
%             end
%             saveas(gcf,[Filename, '.fig'])
%             set(gcf, 'PaperUnits', 'inches');
%             set(gcf, 'PaperPosition', [0 0 imgWidth imgHeight]);
%             print ('-dpng', '-r100', Filename)
%        end
end %for ii=1:32 days per page

L=length(ArrDist(1,1,:));
dum5(1:nTracks,1,1:L)=NaN; dvec5(1:nTracks,1,1:L)=NaN; DistAbs5(1:L)=NaN; DistAbs25(1:L)=NaN;
for ii=1:L
  dum5(1:nTracks-2,1,ii)=abs(ArrDist(1:nTracks-2,1,ii));
  dvec5(ii)=length(find(~isnan(dum5(1:nTracks-2,ii))));
  if dvec5(ii)>50
     DistAbs5(ii)=prctile(abs(ArrDist(1:nTracks-2,1,ii)),5);
     DistAbs25(ii)=prctile(abs(ArrDist(1:nTracks-2,1,ii)),25);
  else
     DistAbs5(ii)=NaN;
     DistAbs25(ii)=NaN; 
  end
end
 

figure('Position',[30 scrsz(4)/20 scrsz(3)/2.3 scrsz(4)/1.15]);
LongA(1:10000)=NaN;
for i=1:length(Long)
    dumLong=Long(i,1:60,ii); dumLat=Lat(i,1:60,ii);
    Vec=find(isnan(dumLong));
    if ~isempty(Vec)
        dumLong(Vec)=[]; dumLat(Vec)=[];
    end
    if length(dumLong)>1
        LongA(i)=interp1(dumLat, dumLong, 29.4);
    end
end
a=hist(LongA, (-15:0.01:-5));


for j=1:floor(nTracks/10)
    i=j*10;
    Vec1=Lat(i,1:60,ii)<29.4;
    Lat(i,Vec1,ii)=29.4;
    Long(i,Vec1,ii)=LongA(i);
    h3=plot(Long(i,1:60,ii),Lat(i,1:60,ii),'-','color',[0.5 0.5 0.5]);
    h3.Color(4) = 0.1;
    hold on
end

load Coast3
for i=1:219
   coastlon=Coast3(i).X; coastlat=Coast3(i).Y;
   plot(coastlon, coastlat, 'k-'); 
   hold on
end

xlim([-15 -5]); ylim([28 38])


if Daynr~=31
    plot([-15 -5], [29.4 29.4], '--','color',rgb('DarkMagenta'),'LineWidth', 1.5);
    plot((-15:0.01:-5), 29.4-smooth(a,100)/(50), '-','color',rgb('DarkMagenta'),'LineWidth', 1.5)
end

plot([-8.7 -5.7 -5.7 -8.7 -8.7],[37.5 37.5 37 37 37.5],'-','color',rgb('Gold'),'LineWidth', 3)

%%determine timing for wind quiver

wegLat=Lat(1:nTracks,1:60,Daynr); 
PassageStep(1:60)=nanmean(wegLat);PassageStep(1)=37.5;
MinLat=PassageStep(1);
for i=1:60    
   if PassageStep(i)>MinLat
      PassageStep(i)=NaN;
   else
       MinLat=PassageStep(i);
   end
end
vec=isnan(PassageStep)==1; PassageStep(vec)=[];

dum(1:60)=nanmean(TIME);
TimeTab=dum(1:length(PassageStep));

Qlat=(30:0.75:37.5); Qlong=(-14.25:0.75:-5.25);LQlat=length(Qlat); LQlong=length(Qlong);
[x,y] = meshgrid(-14.25:0.75:-5.25, 30:0.75:37.5);
Qtime(1:LQlat)=NaN; 
Qwu(1:LQlat,1:LQlong)=NaN;Qwv(1:LQlat,1:LQlong)=NaN;U(1:LQlat,1:LQlong)=NaN;V(1:LQlat,1:LQlong)=NaN;
for i=1:LQlat
    Qtime(i)=interp1(PassageStep, TimeTab, Qlat(i));
    if isnan(Qtime(i)) 
        Qtime(i)=TimeTab(end); 
    end
    for j=1:LQlong
        Qwu(i,j)=interpn(Mlat/100,Mlong/100,MTime, Mu, Qlat(i),Qlong(j),Qtime(i),'linear', -999);
        Qwv(i,j)=interpn(Mlat/100,Mlong/100,MTime, Mv, Qlat(i),Qlong(j),Qtime(i),'linear', -999); 
        DuU=distWB([Qlat(i) Qlong(j)-0.5] ,[Qlat(i) Qlong(j)+0.5]); %km per degree
        DvU=distWB([Qlat(i)-0.5 Qlong(j)] ,[Qlat(i)+0.5 Qlong(j)]); %km per degree
        U(i,j)=Qwu(i,j)/DuU;
        V(i,j)=Qwv(i,j)/DvU;

    end
end

set(gca,'FontSize',16)
dum = gca;
dum.YTick=[28 30  32  34  36  38]; dum.XTick=[-14, -12, -10, -8, -6];
xlabel('Longitude'); ylabel('Latitude'); ylim([28 38]); xlim([-15 -5]);