%% SAPS paper figure: 3x3: SAPS identified at 9 moments. WID with boundary.
% load saved SAPS structure and WID
% Actual Figure 5 in the paper.
close all; clear;
% addpath(genpath('~binzheng/matlab_source'));
% addpath('./subaxis/');
model_type = 'MIX';

%% Directory and runs
rootdir = '/glade/scratch/ldong/'; subdir = 'LTR/';
srun = 'LTR_qq_W17';
hdf_fldr = [pwd,'/'];
% plt_dir0 = ['/glade/work/ldong/Plots/',subdir,srun,'/'];
% plt_dir1 = [plt_dir0,'/SAPS_MIX_TIE/'];
plt_dir  = [pwd,'/'];
% if(~exist(plt_dir1,'dir')); mkdir(plt_dir1); end

%% Start and end time.
dmin = 1; nskip=120;
smin0 = datenum('2013-03-16T18-00-00','yyyy-mm-ddTHH-MM-ss')*1440;
t_start= '2013-03-17T05-00-00'; t_stop = '2013-03-17T21-00-00';
smin = datenum(t_start,'yyyy-mm-ddTHH-MM-ss')*1440;
emin = datenum(t_stop ,'yyyy-mm-ddTHH-MM-ss')*1440;
SW_SM_DAT = load_SW_SM_DAT([hdf_fldr,'SW-SM-DAT']); 
IMF_Bx = SW_SM_DAT(:,7); IMF_By = SW_SM_DAT(:,8); IMF_Bz = SW_SM_DAT(:,9);

%% Load saved SAPS magnitude, latitude, width.
wid0=0.5*1;
load(['SAPS_Structure_variables','_WID0',num2str(wid0*1000,'%4.4d'),'.mat']);%,...
%     'smins','efxnpks','efpklats','efeblats','widnpks','widnpklats','widlatus','widlatls');
load(['saved_Vexb_wi.mat']);%,'xi','yi','WID_smins','WID_Vexb_wis','WID_efxns');

%% Constants and parameters: dipole Br etc.
Rion=1.02;                        % Relative radius of the Earth ionosphere
RE=6.35e6;                        % Radius of the earth [m]
Be=0.30e-4;                       % B field at earth surface [T]
% mix_hdf = dir(fullfile(hdf_fldr,['*mix_',datestr(smin0/1440,'yyyy-mm-ddTHH-MM-ss'),'*']));
% mix_fn = fullfile(hdf_fldr,mix_hdf.name);
% mix_data = getMIXData(mix_fn,model_type);
% xi = mix_data.xi; yi = mix_data.yi; 
zi = sqrt(1-xi.^2-yi.^2);
loni = atan2(yi,xi)*180/pi; colati = acos(zi)*180/pi; lati = 90-colati;
loni1d = loni(:,2); lati1d = lati(2,:);
[nloni,nlati] = size(xi);
dtheta=1.0; mtheta=55;            % Define a regular grid (X,Y) for interpolation
thetac=pi*(-mtheta:dtheta:mtheta)/180;    % On this regular grid, grid points are evenly spaced
[T1, T2]=meshgrid(thetac,thetac);  % which is close to the MIX resolution (dtheta should not be too small)
xdip = Rion.*sin(T1); ydip = Rion.*sin(T2); rdip = Rion; zdip = sqrt(rdip.^2-xdip.^2+ydip.^2);
xydip=sqrt(xdip.^2+ydip.^2); costh = zdip./rdip; cosxr=xdip./xydip; sinyr=ydip./xydip;
B_dipole_r = -2*Be*costh./rdip^3; % [T]
hcells = cell(3,3);

%% Plot
    spc = 0.001; pdn = 0.001; mgn = 0.02; lwid=2; font=16;
    Xlim=[-1,1]*1; Ylim=Xlim; qminmax=[0.1,0.95];
    fig=figure;
    for ii=1:9
        irow = floor((ii-0.5)/3)+1; icol = ii-(irow-1)*3;
        imin = WID_smins(ii); strdate = datestr(imin/1440,'yyyy-mm-ddTHH-MM-ss');
        Vexb_wi = WID_Vexb_wis{ii};
        efxn = WID_efxns{ii};
% Find Equatorial boundary of aurora precipitation
        [efxnpk,iefxnpk] = max(efxn,[],2); % max efxn and lat index
        % [I,J] = ind2sub(siz,IND); linearInd = sub2ind(matrixSize, rowSub, colSub)
        iefpkLine = sub2ind(size(xi),(1:nloni)',iefxnpk); % turn lat index into linear index
        efpklat = lati(iefpkLine);
        ratio = 0.1; iefxneb = zeros(nloni,1)+NaN; % SubAuroral location defined as efxn drops to 0.1 peak.
        for ilon=1:nloni
            iefxnsa = find(efxn(ilon,iefxnpk(ilon):end)>efxnpk(ilon)*ratio); % serach from low to high lat.
            iefxneb(ilon) = iefxnsa(end)+iefxnpk(ilon);
%             iefxnsa = find(efxn(ilon,iefxnpk(ilon):end)<efxnpk(ilon)*ratio);
%             iefxnsa = find(efxn(ilon,iefxnpk(ilon):end)<0.09);
%             iefxneb(ilon) = iefxnsa(1)+iefxnpk(ilon)-1;
        end
        iefebLine = sub2ind(size(xi),(1:nloni)',iefxneb);
        efeblat = lati(iefebLine);
% Find maximum westward ion drift velocity
        widnpk = zeros(nloni,1); iwidnpk = zeros(nloni,1); % peak ion drift and its latitude index.
%         wid0 = 1.0; % only >200m/s can be regarded as WID.
        for ilon=1:nloni
            [widnpk(ilon),iwidnpk(ilon)] = max(Vexb_wi(ilon,iefxneb(ilon):end));
        end % get max WID and index for all lon.
        iwidnpk = iwidnpk+iefxneb-1;
        iwipLine = sub2ind(size(xi),(1:nloni)',iwidnpk);
        widnpklat = lati(iwipLine);
% Find FWHM of WID channel
        widlatu = lati(iwipLine); widlatl = widlatu; lowlim = 45; % 
        for ilon=1:nloni
            if widnpk(ilon)>wid0 % valid WID channel 
                widnpkh = max(0.5*widnpk(ilon),wid0);
                ilat=iwidnpk(ilon);
                % if first point poleward is above 0.5 peak, search polewawrd boundary of WID channel.
                while Vexb_wi(ilon,ilat-1)>widnpkh && efxn(ilon,ilat-1)<=efxnpk(ilon)*ratio
                    ilat = ilat-1;
                end
                widlatu(ilon) = interp1(Vexb_wi(ilon,ilat-1:ilat),lati(ilon,ilat-1:ilat),widnpkh,'linear','extrap');
                if(widlatu(ilon)>efeblat(ilon)); widlatu(ilon)=efeblat(ilon); end
                if(widlatu(ilon)<widnpklat(ilon)); widlatu(ilon)=widnpklat(ilon); end
                ilat=iwidnpk(ilon);
                % if first point equatorward is above 0.5 peak, search equatorwawrd boundary of WID channel.
                while Vexb_wi(ilon,ilat+1)>widnpkh && efxn(ilon,ilat+1)<=efxnpk(ilon)*ratio
                    ilat = ilat+1;
                end
%                 widlatl(ilon) = interp1(Vexb_wi(ilon,ilat:ilat+1),lati(ilon,ilat:ilat+1),widnpkh,'linear','extrap');
                widlatl(ilon) = max(interp1(Vexb_wi(ilon,ilat:ilat+1),lati(ilon,ilat:ilat+1),widnpkh,'linear','extrap'),lowlim);
            end
        end
        iwid = (widnpk>wid0&widlatu<=lati(iefebLine)&(0:360)'<270); % before 12MLT is excluded.
        widnpk(~iwid)=NaN; widnpklat(~iwid)=NaN;
        widlatu(~iwid)=NaN; widlatl(~iwid)=NaN;
        
        hh=subaxis(3,3,ii,'Spacing',spc,'Padding',pdn,'Margin',mgn); hold on; box on; axis equal tight off; hcells{ii}=hh;
        set(hh,'position',[0.0210+(icol-1)*0.29,0.6617-(irow-1)*0.3203,0.3173,0.3173]);
        set(gca,'fontsize',font,'LineWidth',lwid); 
        pcolor(xi,yi,Vexb_wi); shading flat; caxis(quantile(Vexb_wi(:),qminmax)); caxis([-1,1]*3);
        plot(xi(iwipLine(iwid)),yi(iwipLine(iwid)),'r.','Linewidth',lwid);
        plot(cosd(widlatl(iwid)).*cosd(loni1d(iwid)),cosd(widlatl(iwid)).*sind(loni1d(iwid)),'k.','Linewidth',lwid);
        plot(cosd(widlatu(iwid)).*cosd(loni1d(iwid)),cosd(widlatu(iwid)).*sind(loni1d(iwid)),'k.','Linewidth',lwid);
        plot(xi(iefpkLine),yi(iefpkLine),'m','Linewidth',lwid);
        plot(xi(iefebLine),yi(iefebLine),'m:','Linewidth',lwid);
        plotion1('k',[40,50 60 70 80],'sat','MLT',1,font-4); hold on; view(-90,90);
        quiver([-0.8;-0.7],[0.7;0.8],[0.2;0],[0;-0.2],'k','autoscale','off'); 
        quiver(-0.7,0.7,IMF_Bz(imin-smin0-20)/50,-IMF_By(imin-smin0-20)/50,'k','autoscale','off','LineWidth',lwid); 
        text(-0.77,0.58,'B_Y','Fontsize',font-4); text(-0.55,0.86,'B_Z','Fontsize',font-4);
        xlim(Xlim); ylim(Ylim);
        if ii==6
        c1=colorbar; title(c1,'WID [km/s]','Fontsize',font-2);  set(c1,'position',[0.9147    0.3483    0.0206    0.3033]);
        end
        if ii==2
            htt=title('2013-Mar-17'); set(htt,'position',[0.9   -0.0000    0.0000]);
        end
        ht1=text(0.7,0.9,strrep(strdate(end-7:end-3),'-',':'),'Fontsize',font-4);
    end
set(fig,'position',[937    35   888   804]);
% print(fig,[plt_dir,'fig04_WID_evolution','_WID0',num2str(wid0*1000,'%4.4d'),'_r'],'-depsc');
% close(fig);

%% End of code