% 
%  Load relevant best track file, filter for years and months, and apply
%  same spatial filter and vcrit value as used in synthetic event set.
%
%  Copyright WindRiskTech 2012
%  Last modified July 2023 
%--------------------------------------------------------------------------
params
%
load temp vcrit bas shape city_radius clat clong  % Load needed parameters from synthetic track set
xsize=1;
%
bas1=lower(bas);
if strcmp(bas,'AL')
    bas1='at';
end    
file=strcat(trkdir,bas1,'tracks');
%
%  Use Bureau of Meteorology tracks in Australian region
%      (Discontinued in 2013; using IBTrACS for all tracks)
%
%if strcmp(bas,'SH')
%    if max(longstore(:,1)) <= 200 && min(longstore(:,1)) >= 80
%        bas = 'BM';
%        file=strcat(trkdir,'bomtracks');
%    end
%end
%
try
    load(file)
catch err
    msgbox('No best tracks for this basin')
    return 
end   
%
[nnb,mb]=size(vsm);
if strcmp(bas,'BOM')  % Since Bureau uses 10-minute winds, find winds from pressures
    vsm=p_to_v(pcm);
    [utb,vtb,~]=utransbest(latm,longm,monthm,daym,hourm);
    vsm=vsm+0.6.*(1852./3600).*sqrt(utb.^2+vtb.^2); % approximate correction for translation
end   
%
% Reset best track year range if automatic (added June, 2019)
%
if min(bestyears) == 0
    bestbeg=max(min(yeari),min(yearstore));
    bestend=min(max(yeari),max(yearstore));
    if bestend > bestbeg
        bestyears=bestbeg:bestend;
    else
        bestyears=1979:max(yeari);   % Default range of best track years
    end
end  
bestyears=double(bestyears);
%
%  Convert best track wind speeds to knots
%
vsm=vsm.*3600/1852;
%
% If desired, add white noise to IBTrACS wind speeds
%
if strcmpi(add_bt_noise,'y')
    vsm=vsm+windnoise*(0.5-rand(nnb,mb));
end    
%
%  Storm lifetime maximum wind speed
%
vsmmax=max(vsm,[],2);
vsmmax=transpose(vsmmax);
%
%  Condition best track central pressures
%
pcm(pcm==0)=1013;   
%
%  Now apply spatial filters, as with event set
%
passib=zeros(nnb,mb)+1;
passcb=zeros(1,nnb)+1;
passeb=passcb;
if strcmp(shape,'poly')
    [polytype, passib, passeb, xintb, yintb, jintb, kintb, kfracb] = polyfilter(latm,longm);
    %
    % Find maximum ground-relative wind at time track crosses line segments
    %
    [~,xsize]=size(xintb);
    vsegb=zeros(nnb,xsize); % Added 11/28/17
    for nx1=1:xsize % Find intersection with maximum wind
        for i=1:nnb
            %vsegb(i,nx1)=passeb(i).*vsm(i,max(kintb(i,nx1),1)); % Changed 11/28/17
            vsegb(i,nx1)=passeb(i).*(1-kfracb(i,nx1))*vsm(i,max(kintb(i,nx1),1))+ ...  % Changed 7/28/2023
                kfracb(i,nx1).*vsm(i,max(kintb(i,nx1),1)+1);
        end
    end    
    vmaxb=max(vsegb',[],1); %#ok<UDIM>
    if strcmp(polytype,'closed')
        vmaxbcl=max(vsm.*passib,[],2);  % Maximum ground-relative wind speed within closed polygon
        vmaxbcl=transpose(vmaxbcl);        
        vmaxb=max(vmaxb,vmaxbcl);       % For narrow closed polygons, max wind may be on boundary
        clear vmaxbcl
    end    
else  % Case of circular filter, find maximum wind at point of interest
    [passcb, vmaxb, ~] = best_filter(vsm,latm,longm,city_radius,clat,clong);
    if city_radius < 2000 % Exclude cases where city radius encompasses whole basin
        rmb=1540.*max(vsm,1).^(-0.894); % Empirical radius of maximum winds
        rmb=min(rmb,150);
        uinb=zeros(nnb,mb);
        dum=uinb;
        yearbig=repmat(yeari',[1,mb]);
        datem=datenum(double(yearbig),monthm,daym,hourm,0,0); %#ok<DATNM>
        [vsb, dirb, daykb] = pointseriesn(latm,longm,vsm,rmb,dum,dum,datem, ...
                 uinb,uinb,clat,clong,timeres);     
        vsmaxb=max(vsb,[],2);
        vsmaxb=transpose(vsmaxb);
    else
        vmaxb=vsmmax;
    end    
end
[mf,nf]=size(vmaxb);
if mf > nf
    vmaxb=transpose(vmaxb);
end    
%
% Apply time filters
%
passy=ismember(yeari,bestyears);
[ntemp,mtemp]=size(passy);
if ntemp > mtemp 
    passy=transpose(passy);
end    
%
[~,jmin]=min(max((startv-vsm),0),[],2);
passm=zeros(1,nnb);
for i=1:nnb
    passm(i)=ismember(monthm(i,jmin(i)),filtmonths);
end     
%
%  Now find intersection of all filters
%
passt=max(passib,[],2);
[nf,mf]=size(passt);
if nf > mf
    passt=transpose(passt);
end   
passt=max(passt,passeb);
pass=passt.*passcb.*passy.*passm;
if strcmp(shape,'poly') && strcmp(polytype,'open')
    pass=pass.*passeb;
end    
%
%  Maximum velocity filter
%
passv=ismember(round(vmaxb),str2double(vcrit):220);
%
pass=pass.*passv;
%
%  Create filtered variables (Added August, 2021)
%
yearbestt=zeros(1,nnb);
nsib=zeros(1,nnb);
vsmmaxt=zeros(1,nnb);
vmaxbt=zeros(1,nnb);
vsmaxbt=zeros(1,nnb);
if strcmp(shape,'poly')
    vsegbt=zeros(nnb,xsize); % Added 11/28/17
end
monthbestt=zeros(nnb,mb);
daybestt=zeros(nnb,mb);
hourbestt=zeros(nnb,mb);
latbestt=zeros(nnb,mb);
longbestt=zeros(nnb,mb);
vbestt=zeros(nnb,mb);
pbestt=zeros(nnb,mb)+1013;
%
stbest=cell(nnb,1);
[nbestname,~]=size(stname);
ib=0;
for i=1:nnb
    if pass(i) == 1
        ib=ib+1;
        yearbestt(ib)=yeari(i);
        nsib(ib)=nsi(i);
        vsmmaxt(ib)=vsmmax(i);
        vmaxbt(ib)=vmaxb(i);
        if strcmp(shape,'poly')
            vsegbt(ib,:)=vsegb(i,:); % Added 11/28/17
        end
        monthbestt(ib,:)=monthm(i,:);
        daybestt(ib,:)=daym(i,:);
        hourbestt(ib,:)=hourm(i,:);
        latbestt(ib,:)=latm(i,:);
        longbestt(ib,:)=longm(i,:);
        vbestt(ib,:)=vsm(i,:);
        pbestt(ib,:)=pcm(i,:);
        if strcmp(shape,'circ') && city_radius < 2000
            vsmaxbt(ib)=vsmaxb(i);
        end    
        if nbestname > 1
            if nbestname == nnb
                stbest(ib,1)=cellstr(stname(i,:));
            else    
                stbest(ib,1)=cellstr('No Name');
            end    
        else
            stbest(ib,1)=cellstr(stname(:,i));
        end    
    end
end 
if ib < nnb
    vsmmaxt(ib+1:end)=[];
    nsib(ib+1:end)=[];
    vmaxbt(ib+1:end)=[];
    vsmaxbt(ib+1:end)=[];
    stbest(ib+1:end,:)=[];
    if strcmp(shape,'poly')
        vsegbt(ib+1:end,:)=[]; % Added 11/28/17
    end
end
clear vsmmax vmaxb vsmaxb vsegb
vsmmax=vsmmaxt;
vmaxb=vmaxbt;
if strcmp(shape,'poly')
    vsegb=vsegbt; % Added 11/28/17
end
if strcmp(shape,'circ') && city_radius < 2000
    vsmaxb=vsmaxbt;
end    
minbest=zeros(nnb,mb);
secbest=zeros(nnb,mb);
yeartemp=repmat(yearbestt,[mb 1]);
yeartemp=transpose(yeartemp);  
%
% Interpolate to 2 hour time resolution
%
datebestt=datenum(yeartemp, monthbestt, daybestt, hourbestt, minbest, secbest);
vbest=zeros(ib,360);
pbest=zeros(ib,360)+1013;
latbest=zeros(ib,360);
longbest=zeros(ib,360);
datebest=zeros(ib,360);
passib=zeros(ib,360)+1;
passeb=zeros(1,ib)+1;
[~,jmax]=min(abs(latbestt),[],2);
jmax=jmax-1;
jmax=max(jmax,1);
jmax=min(jmax,120);
for i=1:ib
    nj=0;
    for j=1:jmax(i)-1
        %
        % Condition longitudes that cross Greenwich Meridian
        %
        if longbestt(i,1)-longbestt(i,j+1) > 200
            longbestt(i,j+1)=longbestt(i,j+1)+360;
        end    
        for k=1:3
            nj=nj+1;
            w=(k-1)./3;
            vbest(i,nj)=vbestt(i,j).*(1-w)+w.*vbestt(i,j+1);
            pbest(i,nj)=pbestt(i,j).*(1-w)+w.*pbestt(i,j+1);
            latbest(i,nj)=latbestt(i,j).*(1-w)+w.*latbestt(i,j+1);
            longbest(i,nj)=longbestt(i,j).*(1-w)+w.*longbestt(i,j+1);
            datebest(i,nj)=datebestt(i,j).*(1-w)+w.*datebestt(i,j+1);
        end
    end    
    vbest(i,nj+1)=vbestt(i,jmax(i));
    pbest(i,nj+1)=pbestt(i,jmax(i));
    latbest(i,nj+1)=latbestt(i,jmax(i));
    longbest(i,nj+1)=longbestt(i,jmax(i));
    datebest(i,nj+1)=datebestt(i,jmax(i));
end  
if strcmp(shape,'poly')
    [polytype, passib, passeb, xintb, yintb, jintb, kintb, kfracb] = polyfilter(latbest,longbest);
end    
%
%  Sort by local and by liftime maximum wind speeds
%
[vsortb,nsortb]=sortrows(transpose(vsmmax),-1);
vsortb=transpose(vsortb);
nsortb=transpose(nsortb);
%
[vsortpb,nsortpb]=sortrows(transpose(vmaxb),-1);
vsortpb=transpose(vsortpb);
nsortpb=transpose(nsortpb);
if strcmp(shape,'circ') && city_radius < 2000
    [vsortcb,nsortcb]=sortrows(transpose(vsmaxb),-1);
    vsortcb=transpose(vsortcb);
    nsortcb=transpose(nsortcb);
end    
%
clear yearbestt monthbestt daybestt hourbestt minbest secbest datebestt yeartemp ...
    longbestt latbestt vbestt dum utb vtb