function [ nsize ] = prepfilterf(filename, fy2, fm2, filter,varargin )
%
%  Mandatory Inputs:
%
%  filename:     Name of event set (without '.zip')
%  fy2:          Years to use. Can be any MATLAB vector. Set to 0 to use
%                all years. Example:  1980:1985.
%  fm:           Months to use. Set to 1:12 to use all months. 
%  filter:       Type of filter to apply: 'none','circular','basin','or 'poly'. 
%
%  Optional Inputs (in order listed below; use '~' to skip field)
%
%  pf:           Poly file name, without the '.in'. Used if the filter is
%                'poly'. This is an ascii file with 4 columns and any
%                number of rows. Each row contains the longitudes and
%                latitudes of the beginning and end points of each line
%                segment. The line segments may be connected or not, and if
%                the last two columns of the last row are the same as the
%                first two columns of the first row, the filter is taken to
%                be a closed polygon. See 'polyus.in' in the scripts
%                subdirectory for an example. The poly file itself must
%                have the extension '.in'. See User's Guide for information
%                on how to create a poly file using Google Earth.
%  plat:         Latitude of the center of the circle if a circular filter
%                is used. 
%  plong:        Longitude of the center of the circle if a circular filter
%                is used. 
%  prad:         Radius of the circle (km) if a circular filter is ued. 
%  basina:       Ocean basin if the filter type is 'basin'. The choices
%                are: 'al','aus','ep','cp','nio','sp','swio','swp', and
%                'wp', for the North Atlantic, eastern North Pacific,
%                Central North Pacific, North Indian Ocean, South Pacific, 
%                Southwest Indian Ocean, Southwest Pacific, and Western
%                North Pacific, respectively. 
%
%  Output:       nsize is the size (number of tracks) of the event set. 
%
%  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.
%
%  Note that this function clears main event set matrices from base
%  workspace.
%
%  Copyright 2014 WindRiskTech L.L.C.
%  Modified June, 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
%--------------------------------------------------------------------------
%
params    %  Load parameters
filtyears=fy2;
filtmonths=fm2;
%
if nargin > 4
    pf=varargin{1};
    if nargin > 5 && nargin < 9
        plat=varargin{2};
        plong=varargin{3};
        prad=varargin{4};
    end
    if nargin == 9
        basina=varargin{5};
    end          
end    
if strcmpi(filter,'none')
    jfilt=1;
elseif strcmpi(filter,'circular')    
    jfilt=2;
elseif strcmpi(filter,'basin')
    jfilt=3;
    basina=varargin{1};
elseif strcmpi(filter,'poly')
    jfilt=4;
end   
if jfilt == 4
    polyfile=pf;
elseif jfilt == 3
    basininit=basina;
    polyfile=['polysuz' basininit];
elseif jfilt == 2
    poilat=plat;
    poilong=plong;
    if poilong < 0
        poilong=poilong+360;
    end    
    poirad=prad;
end       
%
datea=filename;
storm=[];    
%
if jfilt > 2
    polyfile=['polyfiles/',polyfile '.in'];
end
fs=[char(edir),char(datea),char(storm),char('.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
%
rawfile=dir(fs); 
filedate=rawfile.datenum;
pfiledate=736116;     %   Datenumber corresponding to 1 June 2015, when pressure algorithm was changed in fortran codes
%
%file=[char(datea), char(storm), char('.mat')];
%movefile(char(file),'temp.mat')
if strcmpi(datea(end-2:end),'.nc')
    file=[datea storm];
    convertnc(file)
else    
    file=[datea storm '.mat'];
    movefile(file,'temp.mat','f') 
end   
%
if jfilt > 2 && exist(polyfile,'file') == 0
    errordlg('Poly file does not exist','File Error');
    return
end    
%
if jfilt > 2
    warning('off')
    delete poly.in
    copyfile(polyfile,'poly.in')
    %delete temp.mat
    %delete sorted.mat
    warning('on')
end    
%
load('temp.mat') %#ok<LOAD>
longstore=wrapTo360(longstore);  %#ok<NODEF> % Very, very rarely, a TC will cross the Greenwich line from east to west
save temp.mat longstore -append
if exist('freqyear','var') == 0
    freqyear=freq;   %#ok<NODEF> %  For single year event sets, set freqyear equal to freq
end    
%
if exist('T600store','var')
    [router] = roquick(T600store,rmstore,vstore); %#ok<NODEF>
    save temp router -append
end    
%
if exist('sststore','var')
    save temp sststore -append %#ok<USENS>
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); %#ok<NODEF>
    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
%
[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); %#ok<NODEF>
end    
%
%  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;   %#ok<NODEF>
    clear vms
    save temp pstore -append
end
%
if exist('vsestore','var') == 0 || strcmp(seyewall,'n')
    vsestore=zeros(nn,m);
    rmsestore=zeros(nn,m);
end   
%
%  Calculate translation velocity and wind speed increments
%
if exist('u850store','var')
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore,vstore,u850store,v850store); %#ok<NODEF>
    [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,~]=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); %#ok<NODEF>
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
        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=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<NODEF,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 wil be saved)
%
if strcmp(shape,'poly')
    %if strcmp(polytype,'closed')
        passt=max(passi,[],2)';
        passt=max(passt,passe);
        pass=passt.*passm.*passy; %#ok<NASGU>
    %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; %#ok<NODEF>
vpstoret=vpstore; %#ok<NODEF>
vstoret=vstore;
passit=passi;
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);
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; %#ok<NODEF>
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);
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,:);
        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<AGROW>
            v850store(j,:)=v850storet(i,:);  %#ok<AGROW>
        end    
        if exist('yearstore','var')
            yearstore(j)=yearstoret(i);  %#ok<AGROW>
        end    
        if exist('T600store','var')
            T600store(j,:)=T600storet(i,:);  %#ok<AGROW>
            router(j,:)=routert(i,:);
        end    
        if exist('rhstore','var')
            rhstore(j,:)=rhstoret(i,:);  %#ok<AGROW>
        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,:)=[];
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,:)=[]; %#ok<NASGU>
    v850store(j+1:end,:)=[]; %#ok<NASGU>
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
[nsize,~]=size(vstore);  
%
% 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
%
%  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 vcrit %#ok<USENS>
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 rmsestore vsestore 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
%
evalin( 'base', 'clear *store' )  %  This command clears main arrays from base workspace
disp('prepfilter finished')
end

