%  This script reads in the raw tracks files, runs "utransfull" to determine
%  translation speed, adds increments to the circular wind speed based, and 
%  filters the tracks in space according to the contents of poly files, and
%  by month according to the values contained
%  in the vector filtmonths specified in the params.in file.
%  It orders the tracks by maximum lifetime wind and also
%  by maximum wind at a point of interest, and stores the results. 
%
%  Copyright 2012 WindRiskTech L.L.C.
%  Modified January, 2019
%  Modified January, 2020 to include outer radius calculation
%  Modified September 2020 to write out filtyears
%  Modified December, 2020 to accept netcdf files as input
%--------------------------------------------------------------------------
clear u850store v850store T600store yearstore poilat poilong poirad vsestore rmsestore
clear dir
polyfile='polyfiles/poly.in';   % Default poly file
%
params    %  Load parameters
%
% Issue warning if best tracks being filtered by month
%
aa=ismember(1:12,filtmonths);
amin=min(aa);
if amin == 0
    Mhead=[cellstr('WARNING! Best track months will be filtered according to values in filtmonths in the file params.m.'); ... 
        cellstr('You are advised to reset these to 1:12')];
    warndlg(Mhead,'!!WARNING!!')
end    
%
% Input event set and specify filters
%
Mhead=[cellstr('Note: Years and months will be filtered according to values in filtmonths and filtyears in the file params.m.'); ... 
    cellstr('Spatial filters:')];
UIControl_FontSize_bak = get(0, 'DefaultUIControlFontSize');
set(0, 'DefaultUIControlFontSize', 12);
%{
prompt = {'Enter file name:'};
dlg_title = 'File Name';
num_lines = 1;
answer = inputdlg(prompt,dlg_title,num_lines);
if isempty(answer) || isempty(answer{1})
    return
end    
datea=answer{1};
%}
datea=input('Enter file name ','s');
if datea(1) == '2'
    storm=input('Input storm name (e.g. AL11)','s');
else
    storm=[];
end    
%
jfilt=menu(Mhead, 'None','Circular','By basin','By poly file');
set(0, 'DefaultUIControlFontSize', UIControl_FontSize_bak);
%
if jfilt == 4
    UIControl_FontSize_bak = get(0, 'DefaultUIControlFontSize');
    set(0, 'DefaultUIControlFontSize', 12);
    prompt = {'Input name of poly file (e.g. polyne):'};
    dlg_title = 'Poly File Name';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return
    end    
    polyfile=answer{1};
    set(0, 'DefaultUIControlFontSize', UIControl_FontSize_bak);
elseif jfilt == 3
    set(0, 'DefaultUIControlFontSize', 12);
    prompt = {'Input basin initials (al,ep,wp,nio,swio,aus,or swp):'};
    dlg_title = 'Basin Name';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return 
    end    
    basininit=answer{1};
    set(0, 'DefaultUIControlFontSize', UIControl_FontSize_bak);
    polyfile=['polysuz' basininit];
elseif jfilt == 2
    set(0, 'DefaultUIControlFontSize', 12);
    prompt = {'POI Latitude:','POI Longtitude:','Circle radius (km):'};
    dlg_title = 'POI data';
    num_lines = 1;
    answer = inputdlg(prompt,dlg_title,num_lines);
    if isempty(answer) || isempty(answer{1})
        return 
    end    
    poilat=str2double(answer{1});
    poilong=str2double(answer{2});
    if poilong < 0
        poilong=poilong+360;
    end    
    poirad=str2double(answer{3});
    set(0, 'DefaultUIControlFontSize', UIControl_FontSize_bak);
end    
%
if jfilt > 2
    polyfile=char(['polyfiles/',polyfile,'.in']);
end
fs=[edir,datea,storm,'.zip'];
try
    if exist('7z.exe','file')
        comm=['7z e -aoa -bso0 ',fs];
        system(comm);
    else
        unzip(fs)
    end     
catch err
    errordlg('Unzip error, file may not exist');
    return
end      
%
%  Determine date file was created for use in correcting surface pressure
%
clear dir
rawfile=dir(fs);
try
    filedate=rawfile.datenum;
catch
    disp('Cannot find file')
    return
end    
pfiledate=736116;     %   Datenumber corresponding to 1 June 2015, when pressure algorithm was changed in fortran codes
%
if jfilt > 2 && exist(polyfile,'file') == 0
    errordlg('Poly file does not exist','File Error');
    return
end    
%
clear temp
if jfilt > 2
    warning('off')
    delete poly.in
    copyfile(polyfile,'poly.in','f')    
end    
if strcmpi(datea(end-2:end),'.nc')
    file=[datea storm];
    convertnc(file)
else    
    file=[datea storm '.mat'];
    movefile(file,'temp.mat','f') 
end   
%
clear vcrit
load('temp.mat')
if exist('vcrit','var') == 0 || iscellstr(vcrit) == 0   %#ok<ISCLSTR>
    vcrit='40';  % (default value used for old event sets)
    save temp vcrit -append
end 
longstore=wrapTo360(longstore);  % Very, very rarely, a TC will cross the Greenwich line from east to west
save temp.mat longstore -append
%
[nn,m]=size(vstore);
%
% Following block of code added 10/10/2025 to account for re-scaling of
% rmstore in FORTRAN code.
%
FileInfo = dir('./stats.txt');
TimeStamp = FileInfo.datenum;
if TimeStamp < datenum([2025,10,9,10,0,0]) %#ok<DATNM> 
    rmstore=0.45*rmstore.*sin(pi.*max(latstore,20)./180)./sin(pi*20/180);
end    
%
if exist('freqyear','var') == 0
    freqyear=freq;   %  For single period event sets, set freqyear equal to freq
    save temp freqyear -append
end    
%
if exist('yearstore','var') == 0
    yearstore=zeros(1,nn)+1;    %  For single period event sets, set yearstore = 1
    save temp yearstore -append
end
%
if exist('T600store','var')
    [router] = roquick(T600store,rmstore,vstore);
    save temp router -append
end    
%
if exist('sststore','var')
    save temp sststore -append
end  
%
if jfilt > 2
    shape='poly';
end    
if jfilt == 2
    shape='circ';
    clat=poilat;
    clong=poilong;
    city_radius=poirad;
end    
if jfilt == 3
    if strcmp(basininit,'al') || strcmp(basininit,'ep') || strcmp(basininit,'wp')
        bas=upper(basininit);
    elseif strcmp(basininit,'nio')
        bas='IO';
    else
        bas='SH';
    end    
end
if exist('yearstore','var')
    if filtyears(1) == 0
        filtyears=min(yearstore):1:max(yearstore);
    end    
end
%
% Issue warning of no overlap between years in event set and filtyears in
% params.m
%
if exist('yearstore','var')
    if max(yearstore) < min(filtyears) || min(yearstore) > max(filtyears)
        errordlg('Error: no overlap between years in event set and filtyears in params.m')
        return
    end    
end
%
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
%
%  Correct central surface pressures for pre June 1 2015 event sets
%
if filedate < pfiledate
    vms=vstore.*1852/3600;  %  Convert wind speed to m/s
    pstore=7+pstore.*exp(0.8.*vms.^2./(287*300))-0.6.*vms;  
    clear vms
    save temp pstore -append
end
%
if exist('vsestore','var') == 0 || strcmp(seyewall,'n')
    vsestore=zeros(nn,m);
    rmsestore=zeros(nn,m);
    save temp vsestore rmsestore -append
end   
%
%  Super-gradient wind pre-processor
%
vtemp=vstore;
vupper=105;
vlower=65;
mag=0.25;
%for i=1:nn,
%    for j=1:jmax(i);
%        if vstore(i,j) >= vlower && vstore(i,j) < vupper
%            coef=4.*(vstore(i,j)-vlower).*(vupper-vstore(i,j))./(vupper-vlower)^2;
%            vtemp(i,j)=(1+mag*coef).*vstore(i,j);
%        end
%    end
%end
vstore=vtemp;
clear vtemp
%
%  Calculate translation velocity and wind speed increments
%
if exist('u850store','var')
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore,vstore,u850store,v850store);
    [utse,vtse,uincse,vincse ] = utransfull(latstore,longstore,vsestore,u850store,v850store);
else
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore);
    utse=ut;
    vtse=vt;
    uincse=uinc;
    vincse=vinc;
end
%
vnet=vstore+sqrt(uinc.^2+vinc.^2); % Maximum ground-relative wind speed
[vnetmax,ivmax]=max(vnet,[],2);    % Storm lifetime maximum wind
vnetmax=transpose(vnetmax);
%
% Create mask for secondary eyewall events and find maximum ground relative
% wind for them
%
semask=min(vsestore,1);
vnetse=semask.*(vsestore+sqrt(uincse.^2+vincse.^2));
[vnetmaxse,ivsemax]=max(vnetse,[],2);    % Storm lifetime maximum wind
vnetmaxse=transpose(vnetmaxse);
%
%  Determine genesis points
%
[~,jmin]=min(max((startv-vnet),0),[],2);
%
% Apply time filters
%
passy=1;
if exist('yearstore','var')
    passy=ismember(yearstore,filtyears);
    [nf,mf]=size(passy);
    if nf > mf
        passy=transpose(passy);
    end    
end
passm=zeros(1,nn);
for i=1:nn
    passm(i)=ismember(monthstore(i,jmin(i)),filtmonths);
end        
passi=zeros(nn,m);
%
%  Apply space filters
%
if strcmp(shape,'poly')
    [polytype, passi, passe, xint, yint, jint, kint, kfrac] = polyfilter(latstore,longstore);
    if strcmp(polytype,'closed')
        poly=load('poly.in'); 
        x=poly(:,1); 
        y=poly(:,2);
        A=polyarea(x',y');
        if A < poly_area
            polytype='open';
        end
    end            
    %
    % Find maximum ground-relative wind at time track crosses line segements
    %
    [~,xsize]=size(xint);
    vseg=zeros(nn,xsize); % Added 11/28/17
    for nx=1:xsize % Find intersection with maximum wind
        for i=1:nn
            %vseg(i,nx)=passe(i).*vnet(i,max(kint(i,nx),1)); % Changed 11/28/17
            vseg(i,nx)=passe(i).*(1-kfrac(i,nx))*vnet(i,max(kint(i,nx),1))+ ...  % Changed 7/28/2023
                kfrac(i,nx).*vnet(i,max(kint(i,nx),1)+1);

        end
    end   
    [vmax,iimax]=max(vseg,[],2);
    kmax2=zeros(1,nn);
    for i=1:nn
        kmax2(i)=kint(i,iimax(i));
    end    
    vmax=vmax';  
    if strcmp(polytype,'closed')
        [vmaxcl,kmax2]=max(vnet.*passi,[],2);  % Maximum ground-relative wind speed within closed polygon
        vmaxcl=transpose(vmaxcl);        
        vmax=max(vmaxcl,vmax);         % For narrow closed polygons, max wind may be on boundary
        clear vmaxcl
    end    
else  % Case of circular filter, find maximum wind at point of interest
    if city_radius < 2000 % Exclude cases where city radius encompasses whole basin
        [passc, vmax, kmax2] = best_filter(vnet,latstore,longstore,city_radius,clat,clong);
        yearbig=repmat(yearstore',[1,m]);
        datestore=datenum(double(yearbig),monthstore,daystore,hourstore,0,0); %#ok<DATNM>
        [vs, dir1, dayk] = pointseriesn(latstore,longstore,vstore,rmstore.*magfac, ...
                 vsestore,rmsestore.*magfac,datestore,uinc,vinc, ...
                 clat,clong,timeres);     
        vsmax=max(vs,[],2);
        vsmax=transpose(vsmax);
    else
        vmax=vnetmax;
        kmax2=ivmax;
    end    
end
[mf,nf]=size(vmax);
if mf > nf
    vmax=transpose(vmax);
end    
%
%  Create reduced event set (pass=1 indicates events that will be saved)
%
if strcmp(shape,'poly')
    %if strcmp(polytype,'closed')
        passt=max(passi,[],2)';
        passt=max(passt,passe);
        pass=passt.*passm.*passy;
    %else
    %    pass=passe.*passm.*passy;
    %end
elseif city_radius < 2000
    pass=passc.*passm.*passy;
else
    pass=passm.*passy;
end    
%
% Issue error statement if no events passed filter
%
if max(pass) == 0
    errordlg('Error: No events in the specified event set passed through the specified filter. Please try again')
    return
end    
%
freq=freq.*sum(pass)./nn;  % (In cases with multiple years, this will be revised later in script)
uinct=uinc;
vinct=vinc;
utt=ut;
vtt=vt;
vnett=vnet;
vmaxt=vmax;
vnetmaxt=vnetmax;
uinctse=uincse;
vinctse=vincse;
uttse=utse;
vttse=vtse;
vnettse=vnetse;
vnetmaxtse=vnetmaxse;
daystoret=daystore;
hourstoret=hourstore;
latstoret=latstore;
longstoret=longstore;
monthstoret=monthstore;
pstoret=pstore;
rmstoret=rmstore;
shearstoret=shearstore;
vpstoret=vpstore;
vstoret=vstore;
passit=passi;
kmax2t=kmax2;
if strcmp(shape,'poly')
    passet=passe;
    xintt=xint;
    yintt=yint;
    jintt=jint;
    kintt=kint;
    kfract=kfrac;
elseif city_radius < 2000
    vsmaxt=vsmax;   
    vst=vs;
    dirt=dir1;
    daykt=dayk;
end    
if exist('yearstore','var')
    yearstoret=yearstore;
    yearminorig=min(yearstore);
    yearmaxorig=max(yearstore);
end    
if exist('u850store','var')
    u850storet=u850store;
    v850storet=v850store;
end    
if exist('T600store','var')
    T600storet=T600store;
    routert=router;
end    
if exist('rhstore','var')
    rhstoret=rhstore;
end    
vsestoret=vsestore;
rmsestoret=rmsestore;
j=0;
vmax=zeros(1,nn);
vnetmax=zeros(1,nn);
vnet=zeros(nn,m);
uinc=vnet;
vinc=vnet;
ut=vnet;
vt=vnet;
vnetmaxse=zeros(1,nn);
vnetse=zeros(nn,m);
uincse=vnetse;
vincse=vnetse;
utse=vnetse;
vtse=vnetse;
daystore=zeros(nn,m);
hourstore=daystore;
monthstore=daystore;
latstore=daystore;
longstore=daystore;
rmstore=daystore;
vstore=daystore;
vpstore=daystore;
pstore=daystore;
shearstore=daystore;
passi=zeros(nn,m);
kmax2=zeros(1,nn);
if strcmp(shape,'poly')
    passe=zeros(1,nn);
    [~,xsize]=size(xintt);
    xint=zeros(nn,xsize);
    yint=zeros(nn,xsize);
    jint=zeros(nn,xsize);
    kint=zeros(nn,xsize);
    kfrac=zeros(nn,xsize);
elseif city_radius < 2000
    vsmax=zeros(1,nn);
    [~, vssize]=size(vst);
    vs=zeros(nn,vssize);
    dir1=zeros(nn,vssize);
    dayk=zeros(nn,vssize);
end  
if exist('u850store','var')
    u850store=daystore;
    v850store=daystore;
end
if exist('T600store','var')
    T600store=daystore;
    router=daystore;
end    
if exist('rhstore','var')
    rhstore=daystore;
end 
vsestore=daystore;
rmsestore=daystore;
if exist('yearstore','var')
    yearstore=zeros(1,nn);
end    

for i=1:nn
    if pass(i) == 1
        j=j+1;
        vnet(j,:)=vnett(i,:);
        vnetmax(j)=vnetmaxt(i);
        vmax(j)=vmaxt(i);
        uinc(j,:)=uinct(i,:);
        vinc(j,:)=vinct(i,:);
        ut(j,:)=utt(i,:);
        vt(j,:)=vtt(i,:);
        vnetse(j,:)=vnettse(i,:);
        vnetmaxse(j)=vnetmaxtse(i);
        uincse(j,:)=uinctse(i,:);
        vincse(j,:)=vinctse(i,:);
        utse(j,:)=uttse(i,:);
        vtse(j,:)=vttse(i,:);
        daystore(j,:)=daystoret(i,:);
        hourstore(j,:)=hourstoret(i,:);
        latstore(j,:)=latstoret(i,:);
        longstore(j,:)=longstoret(i,:);
        monthstore(j,:)=monthstoret(i,:);
        pstore(j,:)=pstoret(i,:);
        rmstore(j,:)=rmstoret(i,:);
        shearstore(j,:)=shearstoret(i,:);
        vpstore(j,:)=vpstoret(i,:);
        vstore(j,:)=vstoret(i,:);
        passi(j,:)=passit(i,:);
        kmax2(j)=kmax2t(i);
        if strcmp(shape,'poly')
            passe(j)=passet(i);
            xint(j,:)=xintt(i,:);
            yint(j,:)=yintt(i,:);
            jint(j,:)=jintt(i,:);
            kint(j,:)=kintt(i,:);
            kfrac(j,:)=kfract(i,:);
        elseif city_radius < 2000
            vsmax(j)=vsmaxt(i); 
            vs(j,:)=vst(i,:);
            dir1(j,:)=dirt(i,:);
            dayk(j,:)=daykt(i,:);
        end
        if exist('u850store','var')
            u850store(j,:)=u850storet(i,:); %#ok<SAGROW>
            v850store(j,:)=v850storet(i,:); %#ok<SAGROW>
        end    
        if exist('yearstore','var')
            yearstore(j)=yearstoret(i); 
        end    
        if exist('T600store','var')
            T600store(j,:)=T600storet(i,:); %#ok<SAGROW>
            router(j,:)=routert(i,:);
        end    
        if exist('rhstore','var')
            rhstore(j,:)=rhstoret(i,:); %#ok<SAGROW>
        end    
        vsestore(j,:)=vsestoret(i,:); 
        rmsestore(j,:)=rmsestoret(i,:); 
    end
end
vnet(j+1:end,:)=[];
vmax(j+1:end)=[];
vnetmax(j+1:end)=[];
uinc(j+1:end,:)=[];
vinc(j+1:end,:)=[];
ut(j+1:end,:)=[];
vt(j+1:end,:)=[];
vnetse(j+1:end,:)=[];
vnetmaxse(j+1:end)=[];
uincse(j+1:end,:)=[];
vincse(j+1:end,:)=[];
utse(j+1:end,:)=[];
vtse(j+1:end,:)=[];
daystore(j+1:end,:)=[];
hourstore(j+1:end,:)=[];
latstore(j+1:end,:)=[];
longstore(j+1:end,:)=[];
monthstore(j+1:end,:)=[];
pstore(j+1:end,:)=[];
rmstore(j+1:end,:)=[];
shearstore(j+1:end,:)=[];
vpstore(j+1:end,:)=[];
vstore(j+1:end,:)=[];
passi(j+1:end,:)=[];
kmax2(j+1:end)=[];
if strcmp(shape,'poly')
    passe(j+1:end)=[];
    xint(j+1:end,:)=[];
    yint(j+1:end,:)=[];
    jint(j+1:end,:)=[];
    kint(j+1:end,:)=[];
    kfrac(j+1:end,:)=[];
elseif city_radius < 2000
    vsmax(j+1:end)=[];
    vs(j+1:end,:)=[];
    dir1(j+1:end,:)=[];
    dayk(j+1:end,:)=[];
end    
if exist('u850store','var')
    u850store(j+1:end,:)=[];
    v850store(j+1:end,:)=[];
end    
if exist('T600store','var')
    T600store(j+1:end,:)=[];
    router(j+1:end,:)=[];
end  
if exist('rhstore','var')
    rhstore(j+1:end,:)=[];
end    
vsestore(j+1:end,:)=[];
rmsestore(j+1:end,:)=[];
if exist('yearstore','var')
    yearstore(j+1:end)=[];
end    
%
if exist('yearstore','var')
    yearmin=min(nonzeros(yearstore));
    yearmax=max(yearstore);
    q=max(size(filtyears));
    freqyearnew=zeros(1,q);
    iyn=0;
    for iy=1:q
        if filtyears(iy) >= yearmin && filtyears(iy) <= yearmax
            iyn=iyn+1;
            iyi=filtyears(iy);
            iyiorig=iyi-yearminorig+1;
            iyes=0;
            ino=0;
            for i=1:nn
                if yearstoret(i) == iyi
                    if  pass(i) == 1
                        iyes=iyes+1;
                    else
                        ino=ino+1;
                    end    
                end
            end    
            freqyearnew(iyn)=freqyear(iyiorig)*iyes/(iyes+ino);
        end    
    end 
    if iyn < q
        freqyearnew(iyn+1:end)=[];
    end    
    clear freqyear
    freqyear=freqyearnew;
    freq=mean(freqyear);
end
%
% Clear all temporary arrays
%
clear vmaxt vnett vnetmaxt uinct vinct utt vtt latstoret longstoret yearstoret ...
    monthstoret daystoret hourstoret vstoret vpstoret pstoret rmstoret shearstoret ...
    vsmaxt xintt yintt jintt kintt kfract u850storet v850storet T600storet ...
    passit passet vst dirt daykt vsestoret rmsestoret rhstoret vnettse ...
    vnetmaxset uinctse vinctse uttse vttse kmaxt
%
%  Sort by local and by liftime maximum wind speeds
%
[vsort,nsort]=sortrows(transpose(vnetmax),-1);
vsort=transpose(vsort);
nsort=transpose(nsort);
%
[vsortp,nsortp]=sortrows(transpose(vmax),-1);
vsortp=transpose(vsortp);
nsortp=transpose(nsortp);
%
% Save to files
%
file='sorted';
%
save(file,'vsort','nsort','vsortp','nsortp','freq',...
    'storm','datea')
if strcmp(shape,'circ') && city_radius < 2000
    [vsortc,nsortc]=sortrows(transpose(vsmax),-1);
    vsortc=transpose(vsortc);
    nsortc=transpose(nsortc);
    save sorted vsortc nsortc -append
end    
eset=datea;
if strcmpi(datea(end-2:end),'.nc')
    eset=datea(1:end-3);
end
%
save temp uinc vinc ut vt vnet vmax vnetmax bas  ...
    daystore freq gmeth hourstore latstore longstore monthstore ...
         pstore rmstore shape shearstore vpstore vstore passi kmax2 eset ...
         filtyears vsestore rmsestore vcrit      
if exist('clat','var')
    save temp city_radius clat clong -append
end    
if exist('yearstore','var')
    save temp.mat yearstore freqyear -append
end   
if exist('T600store','var')
    save temp.mat T600store router -append
end         
if exist('rhstore','var')
    save temp.mat rhstore -append
end    
save temp uincse vincse utse vtse vnetse vnetmaxse -append
if strcmp(shape,'poly')
    save temp polytype passe xint yint jint kint kfrac -append
elseif city_radius < 2000
    save temp vs dir1 dayk vsmax -append
end    
%
% Create list of protected variables that will not be cleared by individual
% scripts
%
protected_variables=who('-file', 'temp.mat'); % All variables in 'temp'
%
% Additional protected variables:
%
additional_variables={'protected_variables','x','x2','x3','x4','y','y2','y3', ...
    'y4','z','z2','gplot','tmr','maptoolpresent','xmin*','xmax','ymin','ymax', ...
    'imin','imax','kmin','kmax2', 'histsize', 'sumtotal', 'storm', 'datea', ...
    'fr','nx','ny','nt','rplot','ax','ay','aspect','nplot','hButton'};
%
protected_variables=[protected_variables;additional_variables'];
save temp protected_variables -append
%
disp('prepfilter finished')