%
%  This script reads in the raw tracks files, runs "utransfull" to determinePI
%  translation speed, adds increments to the circular wind speed based,
%  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 November, 2017 (Introduced vseg)
%  Modified January, 2019 to store index of location of maximum winds
%  Modified May, 2019 to use MATLAB's delete and movefile commands
%  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 rhstore rmsestore vsestore freqyear yearstore vs kmax2
clear dir sststore
%
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 name
%
datea=input('Enter file name ','s');
%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};
%
if datea(1) == '2'
    storm=input('Input storm name (e.g. AL11)','s');
else
    storm=[];
end    
%
% Copy event set files to temporary files
%
if exist('poly,in','file')
    delete poly.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 
    errordlg('Unzip error, file may not exist');
    return
end    
%
%  Determine date file was created for use in correcting surface pressure
%
rawfile=dir(fs);
try
    filedate=rawfile.datenum;
catch
    disp('File not found, quitting')
    return
end    
pfiledate=736116;     %   Datenumber corresponding to 1 June 2015, when pressure algorithm was changed in fortran codes
%
if strcmpi(datea(end-2:end),'.nc')
    file=[datea storm];
    convertnc(file)
else    
    file=[datea storm '.mat'];
    movefile(file,'temp.mat','f') 
end    
%
% Load event sets
%
shape='circ';  % (default value used for old event sets)
clear vcrit
load('temp.mat')
if exist('vcrit','var') == 0 || iscell(vcrit) == 0 
    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
%
if exist('freqyear','var') == 0
    freqyear=freq;   %  For single period event sets, set freqyear equal to freq
    save temp freqyear -append
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);
end    
%
% This section of code added on September 27 2020 to deal with tracks over
% Himalayas. This section should not be needed for tracks prior to
% implementation of CD or for after Sept. 27 2020 when FORTRAN code
% corrected
%
[~,jmaxv]=min(vstore,[],2);          % Find when circular wind speeds drop to zero
[~,jmaxlat]=min(abs(latstore),[],2);
%
if sum(jmaxv) < sum(jmaxlat) 
    disp('Himalayan correction implemented')
    for i=1:nn
        latstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        daystore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        hourstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        longstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        monthstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        pstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        rhstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        rmsestore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        rmstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        shearstore(i,jmaxv(i):end)=0;  %#ok<SAGROW>
        T600store(i,jmaxv(i):end)=0;  %#ok<SAGROW>
        u850store(i,jmaxv(i):end)=0; %#ok<SAGROW>
        v850store(i,jmaxv(i):end)=0; %#ok<SAGROW>
        vsestore(i,jmaxv(i):end)=0; %#ok<SAGROW>
        vpstore(i,jmaxv(i):end)=0; %#ok<SAGROW>
    end    
    save temp latstore daystore hourstore longstore monthstore pstore rhstore -append
    save temp rmsestore rmstore shearstore T600store u850store v850store vsestore vpstore -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('sststore','var')
    save temp sststore -append
end    
%
% Define filtyears 
%
filtyears=min(yearstore):max(yearstore); % Use prepfilter if filtering of years desired
%
if exist('T600store','var')
    [router] = roquick(T600store,rmstore,vstore);
    save temp router -append
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;  
    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;
%
%[~,jmax]=min(vstore,[],2);
%jmax=jmax-1;
%
%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);
%
passi=zeros(nn,m);
if strcmp(shape,'poly')
    [polytype, passi, passe, xint, yint, jint, kint, kfrac] = polyfilter(latstore,longstore); %#ok<UNRCH>
    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  
    kmax2=max(kmax2,1);
    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
        [~, 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    
%
%  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 passi kmax2 eset filtyears -append
if strcmp(shape,'poly')
    save temp polytype passe xint yint jint kint kfrac -append %#ok<UNRCH>
elseif city_radius < 2000
    save temp vs dir1 dayk vsmax -append
end    
%
save temp uincse vincse utse vtse vnetse vnetmaxse -append
%
% 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','xb','hButton'};
%
protected_variables=[protected_variables;additional_variables'];
save temp protected_variables -append
%
disp('prep finished')