clear
clc
warning('off')

if 1
    
try
    matlabpool open
end
for yy = 1:90
    try
        close(yy)
    end
end


% physical constants
cnst.m0 = 1.67e-27;  %kg proton mass
cnst.kb = 1.38e-23;  %J/K  Boltzmann constant
cnst.u0 = 4*pi*(1e-7);  %N*UVIS1_trajectory^-2
cnst.const = cnst.u0/2/pi;
cnst.e = 1.602e-19;   %C
cnst.G = 6.67428e-11;  %m3/kg/s2 gravity constant
cnst.m_0 = 5.6846e26;  %kg mass of Saturn
cnst.Rsurf_0 = 56000000;  %m Saturn approximate radius
cnst.m_1 = 1.08022e20;  %kg mass of Enceladus
cnst.Re = 252100;  %m Enceladus radius, assuming perfect sphere
cnst.Rorbit_1 = 237948000;  %m Enceladus orbital radius
cnst.vInit_1 = 2*pi*cnst.Rorbit_1/118386.82;  %m Enceladus orbital velocity
SptlPrcNum = 0;%8;  % these 8 listed first in source crds

% reserve variable space
nmbr = SptlPrcNum + 1;  % 8 Spitale & Porco plus one more
crds = zeros(nmbr,9);
TgrStrpCrds0 = zeros(0,9);

% other perfunctory inputs
sp.SrfSpd = 0;
vp.AngSpd = 0;



%% Model Inputs & Parametrization Resolution
dx = 2000; %m  spacing between tiger stripe points
dxCD = 5000;  %m  intervals for column density integration
xCDMx = 2*100e3;%400e3;%  %m  %max distance to integrate across UVIS trajectory
AnalysisTrj  = [7];  % trajectories to show, UVIS1 (1), UVIS2 (2), E7 (3), E17 (4), E14 (5), E18 (6), E21 (7), E8 (8)
AnalysisTrj2 = [7];  % trajectories to use in Model fit - must be subset of above
k0 = 2;%  %dataset to fit in optimization,  k=1:4 corresponds to masses 2,18,28,44
EnhMx = 20;  % maximum jet enhancement factor over nominal (equal strength) case
dspl = 0;  % update display during optimization
dspint = 1;%200;  % display update interval if dspl = true
repnum = 200;%120;%40;%20;%
JetOptENABLE = 0;  %Perform full optimization?
OSNBanalyze = 1;  %analyze OSNB scans (i.e. E8, E21)??
JetOptFile = 'AvgEnh_E18';%'AvgEnh_UVIS1';%2';%_curtain
LoadPorco = 1;  %loads Carolyn's jets



%%%% SOURCE SETTINGS %%%%
T0 = 130;  %K - inflection point of fit
srcA = 0.0062*0.75*15e6;%(1/0.012)*
srcB = 0;%7;%   %0 is uniform
srcA = srcA*T0^(7-srcB);


Tsrc = 273;  %K  gas temperature - after expanding to collisionless state
TrueAnomFct = [2 0.5 0.5 0.9 0.6 0.5 1 1];% 1];%0*(1:flbylen)+1;% %check last E21 value!!
    TrueAnomFct(:)=1;
SrcRt00 = [1 1 3 0.7];%[0.01 0.01 1 0.01];%
SrcRt00nrm = SrcRt00/sum(SrcRt00);
SrcRt00 = 0.005*1e22*SrcRt00nrm;%(1/0.012)*
Mach00 = [0 2 4 16];%[0 2 5 12];%[0 2 4 12];% 20];%[2 4 12];%  %FOR WATER !!!!

% mSpcs0 = 44*cnst.m0;  %kg  reference species mass
% mSpcs1 = 18*cnst.m0;  %kg  NEW species mass
% MachRatio = sqrt(mSpcs1/mSpcs0);
% Mach00 = MachRatio*Mach00;%Mach00;%


%% Load Tiger Stripes
load TigerStripeCoord.csv  % John Spencer's tiger stripe definition
load TigerStripeTemp.csv  % Carly Howett's tiger stripe temperatures

TgrStrps(1).inds = 1:12;
TgrStrps(2).inds = 14:28;
TgrStrps(3).inds = 30:47;
TgrStrps(4).inds = 49:65;
TgrStrps(5).inds = 67:72;
TgrStrps(6).inds = 74:79;
TgrStrps(7).inds = 81:85;
TgrStrps(8).inds = 87:93;
TgrStrps(9).inds = 95:96;
TgrStrpIntrp = 1;  % makes tiger stripe points evenly spaced in distance
for i = 1:length(TgrStrps)
    % package tiger stripe lats and lons
    TgrStrps(i).lats = TigerStripeCoord(TgrStrps(i).inds,1);
    TgrStrps(i).lons = TigerStripeCoord(TgrStrps(i).inds,2);
    TgrStrps(i).ind = (1:length(TgrStrps(i).inds))';  % indices starting at one
    
    % Tiger Stripe Temperature (K)
    TgrStrps(i).indsT = max(TgrStrps(i).inds,min(TigerStripeTemp(:,1)));
    TgrStrps(i).indsT = min(TgrStrps(i).indsT,max(TigerStripeTemp(:,1)));
    TgrStrps(i).T = interp1q(TigerStripeTemp(:,1), TigerStripeTemp(:,2), TgrStrps(i).indsT');
    %TgrStrps(i).src = 1*2e22*(0*TgrStrps(i).inds + 1)';
    
    % Get Tiger Stripe XYZ
    TgrStrpCrds = TgrStrpCrds0;  % clear coordinates
    sp.lat = TgrStrps(i).lats;
    sp.lon = -TgrStrps(i).lons;
    vp.v0 = 0*sp.lat;
    vp.tlt = 0*sp.lat;
    vp.az = 0*sp.lat;
    TgrStrpCrds = PrtclGenMulti(TgrStrps(i).ind, TgrStrpCrds, cnst.Re, sp, vp);
    TgrStrps(i).crds = TgrStrpCrds(:,1:3);  % Tiger Stripe XYZ
    
    % interpolate tiger stripes
    TgrStrps(i).d = sqrt(sum(((TgrStrps(i).crds(2:end,:) - TgrStrps(i).crds(1:end-1,:)).^2),2));  % distances beteen point - always one less than crds
    TgrStrps(i).dCum = [0; cumsum(TgrStrps(i).d)];  % one less than crds
    TgrStrps(i).len = sum(TgrStrps(i).d);  %m  tiger stripe length
    if TgrStrpIntrp
        dCumMax = (TgrStrps(i).len / dx);
        if dCumMax >= 2
            dCum = (0 : dx: TgrStrps(i).len)';  %m position index along stripe
            dCum = dCum + 0.5*(TgrStrps(i).len - dCum(end));  % center indices
            TgrStrps(i).crds = interp1q(TgrStrps(i).dCum, TgrStrps(i).crds, dCum);
            TgrStrps(i).T = interp1q(TgrStrps(i).dCum, TgrStrps(i).T, dCum);
        end
        
        TgrStrps(i).dCum = dCum;
        TgrStrps(i).d = dCum(2:end) - dCum(1:end-1);
    end
    
    % Tiger Stripe Source Rate -  molecules / m S' Dong et al parameter
    TgrStrps(i).src = srcA*TgrStrps(i).T.^srcB;
    %TgrStrps(i).src = 1*2e22*(0*TgrStrps(i).T + 1)';
    
    TgrStrps(i).num = size(TgrStrps(i).crds,1);
end


%% load INMS data
load Flybys2  %from L1Amulti

% since I saved E14 and E17 in opposite order in L1Amulti - 
% suppose I could also rerun in L1Amulti different order, but whatever
bffr = Flybys(3);
Flybys(3) = Flybys(2);
Flybys(2) = bffr;
clear bffr


%% Make Figure
t0 = -10;%  %s
t1 =  10;%  %s
[JetHnd, hsrf_MDL, trj, dat, Planet] = Enceladus_OSproject3(nmbr, Flybys, AnalysisTrj, cnst, t0, t1);
%[JetHnd, hsrf_MDL, trj, dat, Planet] = Enceladus_OSproject3tst2(nmbr, Flybys, AnalysisTrj, cnst, t0, t1);

%return

%% Set data Cal factors
% make rough estimate for now - get this correct later - BT 050314
OSNBsens = 1/(7.75e15);  % converts 1/m2/s flux to OSNB counts (PER IP!!!)
CO2_H2O = 0.0062;%0.05;%0.0014;%  % CO2 to H2O normalization
flbylen = length(dat);
dat(1).CalFct = (8/6)*9e19/max(dat(1).dat(1).dat);  % quick guess from Hansen et al 2011 - signal to column density
dat(2).CalFct = (15/3.5)*dat(1).CalFct;
dat(3).CalFct = 1;%1/CO2_H2O;%sqrt(293/150)/74.8/((7.11e-10)/1.55)/CO2_H2O/0.031;  % this number times counts per ip gives density in 1/m3
dat(4).CalFct = dat(3).CalFct;
dat(5).CalFct = dat(3).CalFct;
dat(6).CalFct = dat(3).CalFct;
dat(7).CalFct = dat(3).CalFct;
dat(8).CalFct = dat(3).CalFct;
%dat(9).CalFct = dat(3).CalFct;
for yy = 1:flbylen
    for kk = 1:numel(dat(yy).dat)
        dat(yy).dat(kk).dat = dat(yy).CalFct*dat(yy).dat(kk).dat;
    end
end


%% interpolate trajectory & data AND set MDL sigsurf
%dtv = 0.1; %s - for velocity determination only
dt = 0.05;%0.25;%(t1-t0)/8081;% %s
t = (t0:dt:t1)';%(-35:dt:150)';%(1:dt:100)';%
tA = -1e8;  %s  endpoint extension
tB = -tA;  %s  endpoint extension
for i = 1:flbylen
    k0b = min(k0,numel(dat(i).dat));
    vA = (trj(i).trj(2,:) - trj(i).trj(1,:)) / (trj(i).t(2) - trj(i).t(1));
    vB = (trj(i).trj(end,:) - trj(i).trj(end-1,:)) / (trj(i).t(end) - trj(i).t(end-1));
    xA = trj(i).trj(1,:) + vA*tA;  %extend end points to avoid NaNs later
    xB = trj(i).trj(end,:) + vB*tB;  %extend end points to avoid NaNs later
    for k = 1:numel(dat(i).dat)
        dat(i).dat(k).dat0 = dat(i).dat(k).dat;  % record non-interpolated
        dat(i).dat(k).t0 = dat(i).dat(k).t;  % record non-interpolated
            dat(i).dat(k).dat([1,end]) = 0;
        dat(i).dat(k).dat = interp1q([tA;dat(i).dat(k).t;tB], [0;dat(i).dat(k).dat;0], t);
        dat(i).dat(k).dat(isnan(dat(i).dat(k).dat)) = 0;
        dat(i).dat(k).t = t;
    end
    trj(i).trj = interp1q([tA;trj(i).t;tB], [xA;trj(i).trj;xB],t);
    trj(i).ViewTrj = interp1q([tA;trj(i).t;tB], [xA;trj(i).ViewTrj;xB],t);
    %trj(i).trj(isnan(trj(i).trj)) = 1e10;  %m big number
    trj(i).t0 = trj(i).t;  % record non-interpolated
    trj(i).t = t;
    trj(i).inds = (1:size(trj(i).trj,1));
    trj(i).rad = sqrt(sum(trj(i).trj.^2,2));
    [trj(i).radCA, trj(i).indCA] = min(trj(i).rad);
    trj(i).altCA = trj(i).radCA - cnst.Re;
    trj(i).trjNrm = trj(i).trj./(trj(i).rad*[1 1 1]);
        trj(i).vel = (trj(i).trj(3:end,:) - trj(i).trj(1:end-2,:))/(2*dt);
        trj(i).vel = [trj(i).vel(1,:); trj(i).vel; trj(i).vel(end,:)];
%     trj(i).vel = (interp1q([tA;trj(i).t;tB], [xA;trj(i).trj;xB], t+dtv) -...
%         interp1q([tA;trj(i).t;tB], [xA;trj(i).trj;xB], t-dtv))/(2*dtv);
    %trj(i).vel(isnan(trj(i).vel)) = 0;
    trj(i).velMg = sqrt(dot(trj(i).vel,trj(i).vel,2));
    trj(i).d = cumsum(trj(i).velMg*dt);
    for gg = 1:numel(dat(i).dat)  % distanbce scales interpolated to all data sets
        trj(i).d0(gg).d0 = interp1q(trj(i).t,trj(i).d,dat(i).dat(gg).t0);
    end
    %trj(i).d0 = interp1q(trj(i).t,trj(i).d,dat(i).dat(k0b).t0);
    trj(i).velNrm =  trj(i).vel./(trj(i).velMg*[1 1 1]);
    trj(i).velNrm(isnan(trj(i).velNrm)) = 0;
    
    
    % BoreSight direction - Enceladus frame
%     trj(i).BoreSight = 0*trj(i).vel;
%     try
%         trj(i).vINMS=interp1q(trj(i).t0,trj(i).vINMS,trj(i).t);
%         trj(i).BoreSight=interp1q(trj(i).t0,trj(i).BoreSight,trj(i).t);
%         for iii = 1:numel(trj(i).t)
%             trj(i).BoreSight(iii,:) = VctCrdSys(trj(i).vINMS(iii,:), trj(i).vel(iii,:), [-1 0 0]);
%         end
%     end
%     trj(i).BoreSightNrm = trj(i).BoreSight./(sqrt(dot(trj(i).BoreSight,trj(i).BoreSight,2))*[1 1 1]);

    try
        %trj(i).vINMS=interp1q(trj(i).t0,trj(i).vINMS,trj(i).t);
        trj(i).BoreSight = interp1q(trj(i).t0,trj(i).BoreSight,trj(i).t);
        trj(i).BoreSightNrm = trj(i).BoreSight./(sqrt(dot(trj(i).BoreSight,trj(i).BoreSight,2))*[1 1 1]);
        trj(i).Orientation = interp1q(trj(i).t0,trj(i).Orientation,trj(i).t);
    end
    
end
trjlen = size(trj(1).trj,1);



%%  Get UVIS plane points
xCD = -xCDMx:dxCD:xCDMx;  %m  across track displacement distances for column density integration
xCD = xCD + 0.5*(xCDMx - xCD(end));  %center
xCDnum = length(xCD);
for i = 1:flbylen
    % useful trajectory vectors
    r4 = trj(i).trjNrm;%(trj(i).indCA,:);  %surface to peak vector
    r4Mg = sqrt(dot(r4,r4,2));
    r4 = r4./(r4Mg*[1 1 1]);  % make unitary
    if ~trj(i).ViewTrjSET  % if viewing (in UVIS context) from "infinity"
        r1 = trj(i).velNrm;%(trj(i).indCA,:);  %velocity unitary
        r5 = cross(r4,r1,2);  % NOT unitary
    else
        r5 = trj(i).ViewTrj - trj(i).trj;
        %r5 = cross(r4,r6,2);  % NOT unitary
    end
    r5Mg = sqrt(dot(r5,r5,2));
    r5 = r5./(r5Mg*[1 1 1]);  % make unitary - vector perp to trajectory AND surface norm (for UVIS view)
    trj(i).r4 = r4;
    trj(i).r5 = r5;
    
    %trj(i).r5 = (0*trj(i).r5(:,1)+1)*trj(i).r5(end,:);
    
    % get planar points
    trj(i).trjPlane = zeros(trjlen,xCDnum,3);
    for ww = 1:xCDnum
        trj(i).trjPlane(:,ww,:) = trj(i).trj + xCD(ww)*trj(i).r5;%/1000000;
    end
    
    %     if i==1
    %         %line(trj(i).trjPlane(:,:,1),trj(i).trjPlane(:,:,2),trj(i).trjPlane(:,:,3),'Marker','.')
    %         %surface(trj(i).trjPlane(:,:,1),trj(i).trjPlane(:,:,2),trj(i).trjPlane(:,:,3))
    %         line(trj(i).ViewTrj(:,1),trj(i).ViewTrj(:,2),trj(i).ViewTrj(:,3), 'Marker', '.')
    %     end
end
%%% set view to be perpendicular - BT 112816 %%%%%
% i = AnalysisTrj2(1);
% cpsn = trj(i).trj(trj(i).indCA,:)-1e8*trj(i).r5(trj(i).indCA,:);
% set(JetHnd.Parent,'CameraPosition',...
%     [cpsn(1),cpsn(2),cpsn(3)])
% set(JetHnd.Parent,'CameraTarget',trj(i).trj(trj(i).indCA,:))
% set(JetHnd.Parent,'CameraUpVector',trj(i).trj(trj(i).indCA,:))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Other Initializations
MDL = zeros(trjlen,flbylen);
M0 = 5;
N1 = 100;  % resolution of signal surface
N1arr = (1:N1)';
sigsurf_MDL0 = zeros(N1,trjlen,3,flbylen);
JetLen = 0.1e5;  %m Jet display length
for jj = 1:flbylen
    for aa = 1:length(N1arr)
        sigsurf_MDL0(aa,:,:,jj) = permute(trj(jj).trj,[3 1 2]);
    end
end
sigsurf_MDL = sigsurf_MDL0;
trjnum = length(AnalysisTrj);
rat = 0*MDL;
cnt = 1;


%% Get peak properties
CheckPeaks2  % gets peak properties


%%  Mark peak positions along trajectory
% load pkt
figure(1)
MarkPeaks = 0;
if MarkPeaks
    for i = AnalysisTrj
        for ii = 1:length(pkt(i).t)
            pkt(i).hnd(ii) = line(...
                'XData', pkt(i).posdsp(ii,1),...
                'YData', pkt(i).posdsp(ii,2),...
                'ZData', pkt(i).posdsp(ii,3),...
                'Marker', '.',...
                'MarkerSize', 25,...
                'Color', dat(i).clr);
        end
    end
end



%% create Fit diagnostic plot
figure(3)
hplt2 = plot([0],[0],'-','MarkerSize',2);


%% CALCULATE INITIAL PLUME - ON TIGER STRIPE POINTS
MDLPrv = 0*MDL;
MDLbrd = 0*MDLPrv;
AxHnd = get(hsrf_MDL(AnalysisTrj(1)).vl);
AxHnd = AxHnd.Parent;
%Plt = 0;
PorcoSlct = (0*(1:98)+1)';
TgrStrpColor = 1;


% % broad source
% MDLPrv = 0*MDL;
% Mach = 0;
% %SrcRt00 = 2.0*sum(SrcRt00);
% DrawJets = 0;
% Plt = 1;
% DrawFitLine = 0;
% BroadDensity
% MDLbrd = MDLPrv;
% MDLPrv = 0*MDLPrv;

%     Mach = 2;
%     SrcRt0 = SrcRt00(1);
%     DrawJets = 0;
%     Plt = 0;
%     DrawFitLine = 0;
%     EnhanceON = 0;  % whether to alter jet strength
%     LoadPorco = 1;  %loads Carolyn's jets
%     PorcoDensity2
%     SrcRt00(1) = 0.0001*SrcRt00(1);
%     MDLbrd = MDLPrv;
%     MDLPrv = 0*MDLPrv;


% porco jets
DrawFitLine = 0;
EnhanceON = 0;  % whether to alter jet strength
% %%% FRANCIS NIMMO TEST - 112616 %%%%
%     EnhanceON = 1;
%     load PorcoDat%NimmoDat2%NimmoDat
%     AvgEnh = 1*PorcoDat(:,5)+0;%10.^NimmoDat(:,6);%,AnalysisTrj2(1));%
%     clear PorcoDat
%     %AvgEnh = max(AvgEnh,1);
%     AvgEnh = 1.3*AvgEnh/mean(AvgEnh);
%     %AvgEnh(:)=1;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LoadPorco = 1;  %loads Carolyn's jets
DrawJets = 1;%(ijh==1);
TotalSourceArr = 0*SrcRt00;
for ijh = 1:length(SrcRt00)
    Mach = Mach00(ijh);
    SrcRt0 = SrcRt00(ijh);
    SrcRt0nrm = SrcRt00nrm(ijh);
    Plt = (ijh==length(SrcRt00));
    PorcoDensity2
    FitPrms(:,17) = 1;  %enhance factor (see below)
    FitPrmsMach(ijh).FitPrms = FitPrms;  % for different Mach numbers
    FitPrmsMach(ijh).MDL0 = MDL0;
    for ss = 1:size(FitPrms,1)
        [FitPrmsMach(ijh).pkMg(ss,:),BB] = max(MDL0(ss).MDL);  % contribution of jet with equal strength
        FitPrmsMach(ijh).pkt(ss,:) = t(BB)';
    end
    TotalSourceArr(ijh) = TotalSource;
end
% zzz=[(trj(j).d-dCA), MDL(:,j)];
% save S:\SpaceSci\BenT\AAA_TRANSFER\tst zzz
% this source rate DOES NOT include TrueAnomFct - 
% but effect DOES show up in Plot, and
% later optimization COMPENSATES
TotalSourceALL = sum(TotalSourceArr)


%%%% TESTING OLD CODE FOR CAROLYN PLOTS %%%%%
% PorcoEnc = 3;
% AmpSumArr = [];
% jetcnt = 0;
%     ijh = 3;
%     Mach = Mach00(ijh);
%     SrcRt = SrcRt00(ijh);
%     dspCurJet = 0;
% for PorcoSlct = 1:98
%     SrcRt = FitPrmsOpt(:,10);
%     %try
%     PlumeFullDensity2
%     %end
% end
% AmpSumArrSrt = sort(AmpSumArr);
% AmpSumArrSrt = flipud(AmpSumArrSrt);
% AmpSumArrSrt(15)/PorcoEnc
% jetcnt




%%  Evaluate nearby jets
% peak magnitudes
ssarr = (1:size(FitPrms,1))';
MgArr = 0*ssarr*(0*(1:length(pkt))+1);
MgMaxNum = size(FitPrms,1);%  % topmost jets
MgCnd0 = 0*MgArr;
%MgCnd0(1:MgMaxNum,:) = 1;
MgCnd0(end-MgMaxNum+1:end,:) = 1;
for ijh = 1:length(SrcRt00)
    FitPrmsMach(ijh).MgCnd = 0*FitPrmsMach(ijh).pkMg;
    for xx = 1:flbylen
        MgCnd = sortrows([FitPrmsMach(ijh).pkMg(:,xx),ssarr],1);
        MgCnd = sortrows([MgCnd0(:,xx), MgCnd], 3);
        FitPrmsMach(ijh).MgCnd(:,xx) = MgCnd(:,1);
    end
end
%FitPrmsMach(MACH_NUM).pkMg(JET,TRJ)
%FitPrmsMach(MACH_NUM).MgCnd(JET,TRJ)
%FitPrmsMach(MACH_NUM).pkt(JET,TRJ)

% peak times and full condition
ijh = 2;
PkMinDst = 1000*1e3;%100*1e3;%  %m  minimum data peak - to - model distance along trajectory
%cndNumArr = zeros(length(pkt)
for xx = 1:flbylen  %through trjs
    for yy = 1:length(pkt(xx).t)  %through peaks
        PkMdlDst(yy).dst = pkt(xx).velMg(yy)*abs(FitPrmsMach(ijh).pkt(:,xx) - pkt(xx).t(yy));
        PkMdlDst(yy).cndt = (PkMdlDst(yy).dst <= PkMinDst);
        PkMdlDst(yy).cndMg = FitPrmsMach(ijh).MgCnd(:,xx);
        PkMdlDst(yy).cnd = (PkMdlDst(yy).cndt.*PkMdlDst(yy).cndMg == 1);
        PkMdlDst(yy).cndNum = sum(PkMdlDst(yy).cnd);
        cndNumArr(yy) = PkMdlDst(yy).cndNum;
    end
    pkt(xx).cndALL = all(cndNumArr);
    pkt(xx).PkMdlDst = PkMdlDst;
end
%pkt(:).cndALL
%pkt(TRJ).PkMdlDst(PEAKNUM).cnd(JET)




%%
% Go through each peak, check all jets, and enhance ONLY best fit jet
% Influenced by order in which it is done
if JetOptENABLE
    FitPrmsMach0 = FitPrmsMach;
    %enhfctLst = 0*(1:size(FitPrms,1)) + 1;
    DevVlALL = zeros(repnum,1);
    enhfctLstALL = zeros(size(FitPrms,1),repnum);
    for qqq = 1:repnum
        %tic
        
        %% Reset Model
        %MDLPrv = 0*MDL;  % starting plume (zero OR 'diffuse' background)
        
        
        %% Randomly select peaks
        % concatonate peak list
        pklst = zeros(0,2);
        for i = AnalysisTrj2%AnalysisTrj
            %if ~trj(i).dsg  %INMS - only try to fit INMS, use UVIS as a constraint only
                pknm = length(pkt(i).t);
                pklstNew = zeros(pknm,2);
                pklstNew(:,1) = i;
                pklstNew(:,2) = 1:pknm;
                pklst = [pklst; pklstNew];
            %end
        end
        
        % scramble peak list
        pklstRand = zeros(0,size(pklst,2));
        pklstLen = size(pklst,1);
        for i = 1:pklstLen
            rndind = round(rand()*(size(pklst,1)-1))+1;
            pklstRand(i,:) = pklst(rndind,:);
            pklst(rndind,:) = [];
        end
        pklst = pklstRand;
        clear pklstRand
%         figure(20)
%         plot(pklst(:,1),pklst(:,2),'-o')  %of course there are huge number of permutations
        
        % set indices to scrambled peak list
        trjind = pklst(:,1);  % flyby designations
        pkind = pklst(:,2);  % peak designations
        %     trjind = [5, 6, 4];  % flyby designations
        %     pkind = [4, 5, 4];  % peak designations
        %     trjind = [2, 2];  % flyby designations
        %     pkind = [3, 2];  % peak designations
        
        
        
        %% Loop through data peaks
        %%%%%  FitPrms  %%%%%
        %1. opt iterations per peak
        %2. time interval on time scale at CA - used for M fit
        %3. deviation data - model
        %4. swivel angle
        %5. which TS (ind)
        %6. where on TS (ind)
        %7. swivel angle (ind)
        %8. Mach number
        %9. angle of Jet to surface normal
        %10. prfct Jet strength
        %11-13.  Jet unit vector
        
        % eval current signals
        FitPrmsMach = FitPrmsMach0;
        MDL0prv = 0*FitPrmsMach(1).MDL0(1).MDL + MDLbrd;  % baseline profile, NOTE: MODIFIED EVERYTIME A PEAK IS FIT !!
        for ijh = 1:length(SrcRt00)
            for ss = 1:size(FitPrms,1)
                MDL0prv = MDL0prv + FitPrmsMach(ijh).MDL0(ss).MDL;
            end
        end
        
        %FitPrmsMachNew = FitPrmsMach;
        ijharr = 1:length(SrcRt00);
        MnVlTot = 0;
        pklstLen1 = pklstLen;
        enhfctLst = 0*(1:size(FitPrms,1)) + 1;
        for gg = 1:length(trjind)  % scans through peaks - identify best fitting jet (see end)
            
            FitPrmsMachNew = FitPrmsMach;
            PkPosInd = int32(pkt(trjind(gg)).ind(pkind(gg)));  % index of peak position along trajectory
            if MarkPeaks
                MrkrSz = get(pkt(trjind(gg)).hnd(pkind(gg)));
                MrkrSz = MrkrSz.MarkerSize;
                set(pkt(trjind(gg)).hnd(pkind(gg)),'MarkerSize',2*MrkrSz)
            end
            
            
            %% vary Jet STRENGTH and update display
            DevVl = [];
            ssarrCut = ssarr(pkt(trjind(gg)).PkMdlDst(pkind(gg)).cnd)';%ssarr;%
            PkPosIndALL = int32(pkt(trjind(gg)).ind);
            for ss = ssarrCut%ssarr%
                
                % eval current signals
                MDL = 0*FitPrmsMach(1).MDL0(1).MDL;
                for ijh = ijharr
                    MDL = MDL + FitPrmsMach(ijh).MDL0(ss).MDL;
                end
                MDLprv = MDL0prv - MDL;  % get (TEMPORARILY) contribution of baseline profile in other jets
                
                % renorm MDL
                for i = AnalysisTrj
                    k0b = min(k0,numel(dat(i).dat));
                    rat(:,i) = max((dat(i).dat(k0b).dat - MDLPrv(:,i)),0) ./ MDL(:,i);  %enhancement ratio on jet, from its baseline magnitude
                end
                %rat(PkPosInd,trjind(gg))
                enhfct = rat(PkPosInd,trjind(gg));  % max ratio for current trajectory
                    enhfct(enhfct*enhfctLst(ss) > EnhMx) = 1;
                    enhfct(isnan(enhfct)) = 1;
                MDL = enhfct*MDL;
                MDLprv = MDLprv + MDL;  % reverse above, this is MODIFIED SIGNAL
                
                % deviation
                dev = 0*MDL(:,1);
                for i = AnalysisTrj2
                    k0b = min(k0,numel(dat(i).dat));
                    dev = dev + abs(dat(i).dat(k0b).dat - MDLprv(:,i));  % deviation of new jet contribution from that 'required'
                end
                DevVl = [DevVl, sum(dev)];  % accumulate deviations for all jets, to current peak
                
                % display
                if ~mod(cnt,dspint) && dspl
                    %DrawFitLine = 0;
                    PkFitDspl
                end
                cnt = cnt + 1;
                
                % eval current signals
                %MDL = 0*FitPrmsMach(1).MDL0(1).MDL;
                for ijh = ijharr
                    FitPrmsMachNew(ijh).MDL0(ss).MDL = enhfct*FitPrmsMach(ijh).MDL0(ss).MDL;
                end
                
                % record
                % accumulate enhancement factors above baseline for all jets, to current peak
                % RHS has factor of previous enhfctLst(ss), since MDL0prv changed after every peak fit
                enhfctLst(ss) = enhfct*enhfctLst(ss);
            end
            
            
            %% Re-plot Final Result at Best Fit
            
            [MnVl, MnInd] = min(DevVl);
            %cndD = pkt(trjind(gg)).PkMdlDst(pkind(gg)).cnd;%(rand(size(FitPrms,1),1) > 0.1);
            %[MnVl, MnInd] = min(DevVl(cndD==1));
            
            MnInd = ssarrCut(MnInd);
            % eval current signals
            if any(MnInd)
                MDL = 0*FitPrmsMach(1).MDL0(1).MDL;
                for ijh = ijharr
                    MDL0prv = MDL0prv - FitPrmsMach(ijh).MDL0(MnInd).MDL + FitPrmsMachNew(ijh).MDL0(MnInd).MDL;
                    FitPrmsMach(ijh).MDL0(MnInd).MDL = FitPrmsMachNew(ijh).MDL0(MnInd).MDL;
                    FitPrmsMach(ijh).FitPrms(MnInd,17) = enhfctLst(MnInd);  % enhancement factor
                end
                MDL = MDL0prv;
                if gg==length(trjind)
                    %DrawFitLine = 0;
                    %PkFitDspl
                end
            else
                MnVl = 0;
                pklstLen1 = pklstLen1 - 1;
            end
            
            
            %% collect results
            enhfctLst = FitPrmsMach(1).FitPrms(:,17);
            MnVlTot = MnVlTot + MnVl;
            
            
            %%  some display elements
            %gg
%             try
%                 set(JtHndTmp, 'Visible', 'off')
%             end
%             JtHndTmp = JetHnd2(MnInd);
%             set(JtHndTmp, 'Visible', 'on')
            %pause(1)
            if MarkPeaks
                set(pkt(trjind(gg)).hnd(pkind(gg)),'MarkerSize',MrkrSz)
            end
        end
        
        
        %% collect results
        %MnVlTot = MnVlTot/pklstLen1;
        for i = AnalysisTrj2
            k0b = min(k0,numel(dat(i).dat));
            dstThrsh = 150*1e3;  %+/- distance range to consider in deviation calc
            dCA = trj(i).d(trj(i).indCA);  %m  CA position
            cnd = (abs(trj(i).d - dCA) < dstThrsh);
            DevVlALL(qqq) = sum(abs(dat(i).dat(k0b).dat(cnd) - MDL(cnd,i)));%MnVlTot;
        end
        enhfctLstALL(:,qqq) = FitPrmsMach(1).FitPrms(:,17);
        AvgEnh = sum(enhfctLstALL,2)/qqq;  %%% USED BELOW IN AVERAGE RESULT %%%
        
        
        %% plot
        figure(46)
        plot(DevVlALL,'-o')
        %figure(47)
        %plot(AvgEnh,PorcoDat(:,5),'o')
        figure(48)
        ijk2 = 2;
        MeanCntrb = mean(FitPrmsMach(ijk2).pkMg(:,AnalysisTrj2),2);
        MeanCntrb = MeanCntrb / max(MeanCntrb);
        plot(ssarr,FitPrmsMach(1).FitPrms(:,17),'-o',...
            ssarr,MeanCntrb,'-o')
%         plot(AvgEnh,MeanCntrb,'o')
        
        PlottingInstructions2
    end
    %qqq
    
    
    %% Save Average Result
    %save AvgEnh_E17 AvgEnh enhfctLstALL DevVlALL
    save(JetOptFile,'AvgEnh','enhfctLstALL','DevVlALL','repnum')

else
    %load AvgEnh1%_E14%E18
    load(JetOptFile)
        repnum = size(enhfctLstALL,2);  %%OLD
end



%% Solution Distance Matrix, and Categorization
dstMat = zeros(size(enhfctLstALL,2),size(enhfctLstALL,2));
for ii = 1:size(enhfctLstALL,2)
    for jj = 1:size(enhfctLstALL,2)
        dstMat(ii,jj) = sqrt(sum((enhfctLstALL(:,ii) - enhfctLstALL(:,jj)).^2));
    end
end

%categorization routine
cnt2Prv = -1;
cnt2 = 0;
dstStrt = median(dstMat(dstMat>0)) + 2*std(dstMat(dstMat>0));
dstEnd = median(dstMat(dstMat>0)) - 2*std(dstMat(dstMat>0));
dstStrt = max(dstStrt,median(dstMat(dstMat>0))/3);
dstEnd = max(dstEnd,median(dstMat(dstMat>0))/3);
dstStp = std(dstMat(dstMat>0))/3;
dstVl = dstStrt;
while (cnt2 >= cnt2Prv) && (dstVl >= dstEnd)  % loops through distance limiters to find max categories
    cnt2Prv = cnt2;
    TstLst = zeros(size(enhfctLstALL,2),1);
    StrtSum = sum(TstLst);
    cnd00 = (dstMat < dstVl);
    cnd00 = cnd00.*(dstMat>0);
    cnd = 0*cnd00 + 1;
    cnt2 = 0;
    RecImg = [];
    while 1 % loops until all categories found
        cnd0 = cnd.*cnd00;
        TstLst = sum(cnd0,2) > 0;
        TstLst(cumsum(TstLst) > 1) = 0;
        if ~sum(TstLst)
            break
        end
        RecImg0 = 0*TstLst;
        while sum(TstLst)
            RecImg0(TstLst>0) = 1;
            TstLstPrv = (TstLst > 0);
            cnd0(TstLstPrv==1,:) = 0;
            cnd(TstLstPrv==1,:) = 0;
            TstLst = cnd0*TstLst;
        end
        RecImg = [RecImg, RecImg0];
        cnt2 = cnt2+1;
    end
    RecImg = [RecImg,1-sum(RecImg,2)];
    sum(RecImg,2)
    dstVl = dstVl - dstStp
    cnt2
    %pause(0.5)
end
if ~cnt2
    RecImg = 0*(1:repnum)' + 1;
end



%% Plot and Analyze Solution Groupings - Discern Correlaions
figure(55)
TgrStrpColor = 0;
hpltGrp(1) = plot([0],[0],'-','Marker', '.', 'MarkerSize', 25);
hpltGrpPrnt = get(hpltGrp(1));
hpltGrpPrnt = hpltGrpPrnt(1).Parent;
for gg = 2:cnt2+1
    hpltGrp(gg) = copyobj(hpltGrp(1),hpltGrpPrnt);
end
DevAvgLst = 0*(1:cnt2+1);
AvgEnhLst = 0*ssarr*DevAvgLst;
ssarrALL = ssarr*(0*(1:repnum)+1);
TotalSourceArr = 0*DevAvgLst;
for gg = 1:cnt2+1%repnum%
    
    AvgEnh = mean(enhfctLstALL(:,RecImg(:,gg)==1),2);  % plots average of solutions
        AvgEnh(isnan(AvgEnh(:))) = 0;
    AvgEnhLst(:,gg) = AvgEnh;
    DevAvgLst(gg) = mean(DevVlALL(RecImg(:,gg)==1));  % average of deviations of ORIGINAL solutions
        DevAvgLst(isnan(DevAvgLst(:))) = 0;
    %AvgEnh = enhfctLstALL(:,gg);%
    %AvgEnh(:) = 0*ssarr + 1.5;
    %pause(3)
    
    %plot
    GrpEnhVls = enhfctLstALL(:,RecImg(:,gg)==1);
    ssarrVls = ssarrALL(:,RecImg(:,gg)==1);
    PltShft = 0.2*(0*ssarrALL(:,1)+1)*((0:(size(ssarrVls,2)-1))/(size(ssarrVls,2)-1));
    GrpEnhVls = GrpEnhVls + PltShft;
    AA = sortrows([ssarrVls(:), GrpEnhVls(:)]);
    GrpEnhVls = AA(:,2);
    ssarrVls = AA(:,1);
    clr = rand(1,3);%
    set(hpltGrp(gg), 'XData', ssarrVls(:), 'YData',GrpEnhVls(:)+0.1*gg,...
        'Color', clr)
    
    %AvgEnh(:)=1;
    
    %% Re-eval jets with above result
    DrawFitLineEVER = 0;
    MDLPrv = 0*MDL;
    TotalSourceArr = 0*SrcRt00;
    for ijh = 1:length(SrcRt00)
        Mach = Mach00(ijh);
        SrcRt0 = SrcRt00(ijh);
        SrcRt0nrm = SrcRt00nrm(ijh);
        DrawJets = 0;%(ijh==1);
        Plt = (ijh==length(SrcRt00));
        DrawFitLine = DrawFitLineEVER*(ijh==length(SrcRt00));
        EnhanceON = 1;  % whether to alter jet strength
            %LoadPorco = 1;  %loads Carolyn's jets
        PorcoDensity2
        FitPrms(:,17) = 1;  %enhance factor (see below)
        FitPrmsMach(ijh).FitPrms = FitPrms;  % for different Mach numbers
        FitPrmsMach(ijh).MDL0 = MDL0;
        for ss = 1:size(FitPrms,1)
            [FitPrmsMach(ijh).pkMg(ss,:),BB] = max(MDL0(ss).MDL);
            FitPrmsMach(ijh).pkt(ss,:) = t(BB)';
        end
    TotalSourceArr(ijh) = TotalSource;
    end
TotalSourceALL(gg) = sum(TotalSourceArr);
    
end
TotalSourceALL = TotalSourceALL'

% %%%% PLUME FIGURE - Colorize Average of solution groupings %%%%
% Plt = 0;
% DrawJets = 0;
% AvgEnhLst(:,1) = mean(AvgEnhLst,2);
% %     AvgEnhLst(:)=1;
% %     AvgEnhLst(59) = 20;
% gg=1;
% % [~,gg] = min(DevAvgLst);
% %LoadPorco = 1;  %loads Carolyn's jets
% PorcoDensity2
% %%%%

end


%% Extrapolate to new Mach number and compare to data in other mass channels
%load xxx
if 1
% New Species and Mach number
mSpcs0 = 44*cnst.m0;  %kg  reference species mass
mSpcs1 = 18*cnst.m0;  %kg  NEW species mass
Tgas = 150;  %K  jet gas temperature
vthm1 = sqrt(2*cnst.kb*Tgas/mSpcs1);  %m/s  NEW species most probable thermal speed
MachRatio = sqrt(mSpcs1/mSpcs0);
Mach00new = MachRatio*Mach00;%Mach00;%
MixRatio = 1/0.0036/1.3282;%/0.5866/0.2333;%1/0.004/0.9116;  %mixing ratio of new species
    %Mach00new(:) = 16;
    
k0 = 2;  %data select - reset

DrawFitLine = 1;
%EnhanceON = 0;  % whether to alter jet strength
%LoadPorco = 1;  %loads Carolyn's jets
%DrawJets = 0;%(ijh==1);

%AnalysisTrj2 = [5];
TotalSourceALL = 0*(1:numel(DevAvgLst));
MDLprfs = zeros(size(MDL,1),size(MDL,2),numel(DevAvgLst));  % for origin - BT 102716
for gg = 1:numel(DevAvgLst)
AvgEnh = AvgEnhLst(:,gg);
    %AvgEnh(:)=1;
    AvgEnh = mean(AvgEnhLst,2); %%%%
    %AvgEnh = AvgEnhLst(:,1); %%%%
MDLPrv = 0*MDL;
TotalSourceArr = 0*SrcRt00;
MDL0arrALL = zeros([size(MDL0arr),numel(SrcRt00)]);
DrawJets = 0;
for ijh = 1:length(SrcRt00)
    Mach = Mach00new(ijh);
    SrcRt0 = MixRatio*SrcRt00(ijh);
    SrcRt0nrm = SrcRt00nrm(ijh);
    Plt = (ijh==length(SrcRt00));
    DrawJets = (ijh==length(SrcRt00)).*(gg==numel(DevAvgLst));
    PorcoDensity2
    FitPrms(:,17) = 1;  %enhance factor (see below)
    FitPrmsMach(ijh).FitPrms = FitPrms;  % for different Mach numbers
    FitPrmsMach(ijh).MDL0 = MDL0;
    for ss = 1:size(FitPrms,1)
        [FitPrmsMach(ijh).pkMg(ss,:),BB] = max(MDL0(ss).MDL);  % contribution of jet with equal strength
        FitPrmsMach(ijh).pkt(ss,:) = t(BB)';
    end
    TotalSourceArr(ijh) = TotalSource;
    MDL0arrALL(:,:,:,ijh) = MDL0arr;
end
TotalSourceALL(gg) = sum(TotalSourceArr);
% TotalSourceALL(gg) = sum(TotalSourceArr);
% TotalSourceALL = TotalSourceALL'
j = AnalysisTrj2(1);
k0b = min(k0,numel(dat(j).dat));
MixRatio2(gg) = sum(MDL(:,j))/sum(dat(j).dat(k0b).dat);
TotalSourceALL(gg) = TotalSourceALL(gg)/MixRatio2(gg);
SrcIndiv(:,gg) = SrcRts*m/MixRatio2(gg);

%%%  jet source rates, run with mass to 18, k0 to 2
%%% then MixRatio above to 1 over value of MixRatio2, when comparing to H2O
%%% rationale would not apply to UVIS, there MixRatio just affects the fit
JrcInd = 98;
SrcRts(JrcInd)*m/MixRatio2(gg)
%%%

%axis([-250 250 0 0.5e-10])
MDLprfs(:,:,gg) = MDL;
end
TotalSourceALL = TotalSourceALL'
MixRatio2 = MixRatio2'


%%%% for origin %%%%
a1 = AvgEnhLst;
a3 = [TotalSourceALL, MixRatio2, DevAvgLst'];
%a2 = MeanCntrb;
%%%%
%zzz=[trj(j).t, (trj(j).d-dCA)/1000, MDL(:,9)]
trjID = j;
zzz = [trj(j).t, (trj(j).d-dCA), permute(MDLprfs(:,trjID,:),[1 3 2])];
zzzdat = [dat(trjID).dat(k0b).t0, (trj(trjID).d0(k0b).d0 - dCA), dat(trjID).dat(k0b).dat0];
%save S:\SpaceSci\BenT\AAA_TRANSFER\tst a1 a3 zzz zzzdat
end



%%  Heighlight active jets

try
    AvgEnhThrsh = 1.5;
    IndShow = [];
    for ii = 1:length(AvgEnh)
        if AvgEnh(ii) > AvgEnhThrsh
            set(JetHnd2(ii), 'Visible', 'on')
            IndShow = [IndShow, ii];
        else
            set(JetHnd2(ii), 'Visible', 'off')
        end
    end
end


% draw contribution of each individual simulated jet
DrawEachJetPrfs = 0;
if DrawEachJetPrfs
    for rr = 1:size(FitPrms,1)
        if PorcoSlct(rr)%4.6%2.5
            amp = 5.5*(dat(j).sclfct/dat(j).CalFct)*permute(sum(MDL0arrALL(rr,:,j,:),4),[2 3 1]);
            AmpSum = sum(amp);
            if AmpSum >= 0%1e6
                sigsurf_MDL(:,:,3,j) = sigsurf_MDL0(:,:,3,j) - ((N1arr-1)/N1)*amp';
                clr = PorcoClr(rr,:);
                line('XData', sigsurf_MDL(1,:,1,j),...
                    'YData', sigsurf_MDL(1,:,2,j),...
                    'ZData', sigsurf_MDL(1,:,3,j)-amp',...
                    'LineWidth', 3,'Color',clr,...
                    'Tag', num2str(rr),...
                    'Parent',AxHnd)
            end
        end
    end
end

% SnpshtFormat = 'bmp';
% SnpshtFlNm = ['tstC', JetOptFile];%'tst18';
% Fig_image = frame2im(getframe(1));
% imwrite(Fig_image,[SnpshtFlNm,'.',SnpshtFormat],SnpshtFormat);




%% Record Image Slices
%PlumeFullDensityImage102214

%end



%% OSNB Scan Test - BT 040415

if OSNBanalyze
% New surface image
figure(1)
caxis([0 1e4])  % gas velocity range of surface projection color scale
iiiarr = 1:1:numel(trj(i).t);

% Create Enceladus Representation for plot
N = 800;%200;%
Enceladus = zeros(N,N,3);
aa = 1:N;
ang1a = (180-90)*(pi/180);
ang1 = ang1a + (pi-ang1a)*(aa-1)/(N-1);
ang2 = 2*pi*(aa-1)/(N-1);
Planet2(:,:,1) = cnst.Re*sin(ang1)'*cos(ang2);
Planet2(:,:,2) = cnst.Re*sin(ang1)'*sin(ang2);
Planet2(:,:,3) = cnst.Re*cos(ang1)'*(0*ang2+1);
alphaImg0 = 0*Planet2(:,:,1);
%alphaImg(Planet2(:,:,3) < -0.95*cnst.Re) = 1;
% surface projection shape
SrfHnd = surface(1.01*Planet2(:,:,1),1.01*Planet2(:,:,2),1.01*Planet2(:,:,3),...
    'CData',0*Planet2(:,:,1),...
    'LineStyle','none',...
    'FaceColor','interp',...
    'FaceAlpha','interp',...
    'AlphaDataMapping','scaled',...
    'AlphaData',alphaImg0);



% set v_vals and Evals scales
v0 = 10;  %m/s
v1 = 20000;  %m/s
v_fct = 1/50;
v_val = 0:v0*v_fct:v0;
v_new = v_val(end)*(1 + v_fct);
while v_new <= v1
    v_val = [v_val, v_new];
    v_new = v_val(end)*(1 + v_fct);
end
dv = (v_val(3:end)-v_val(1:end-2))/2;
dv(end+1) = v_val(end)-v_val(end-1);
dv = [(v_val(2)-v_val(1)), dv];


% New Species and Mach number
mSpcs0 = 18*cnst.m0;  %kg  reference species mass
mSpcs1 = 18*cnst.m0;%2*cnst.m0;  %kg  NEW species mass
Tgas = 85;%  %K  jet gas temperature
vthm1 = sqrt(2*cnst.kb*Tgas/mSpcs1);  %m/s  NEW species most probable thermal speed
MachRatio = sqrt(mSpcs1/mSpcs0);%sqrt(Tgas/150);%
Mach00new = MachRatio*Mach00;%Mach00;%
MixRatio = (1800/2800)*1*MixRatio;%/0.005;%0.01;%  %mixing ratio of new species
vrt = v_val/vthm1;
vrt2D = (0*PorcoSlct+1)*vrt;
    %SrcRt00(1) = 0;


% Re-eval jets with NEW Mach number
DrawFitLineEVER = 0;
EnhanceON = 0;  % whether to alter jet strength
%LoadPorco = 0;  %loads Carolyn's jets?
DrawJets = 0;
Plt = 1;
DrawFitLine = 0;
MDLPrv = 0*MDL;
MDL0arrALL = zeros([size(MDL0arr),numel(SrcRt00)]);
for ijh = 1:numel(SrcRt00)
    Mach = Mach00new(ijh);
    SrcRt0 = MixRatio*SrcRt00(ijh);
        SrcRt0nrm = SrcRt00nrm(ijh);
    PorcoDensity2
    FitPrms(:,17) = 1;  %enhance factor (see below)
    FitPrmsMach(ijh).FitPrms = FitPrms;  % for different Mach numbers
    FitPrmsMach(ijh).MDL0 = MDL0;
    for ss = 1:size(FitPrms,1)
        [FitPrmsMach(ijh).pkMg(ss,:),BB] = max(MDL0(ss).MDL);
        FitPrmsMach(ijh).pkt(ss,:) = t(BB)';
    end
    MDL0arrALL(:,:,:,ijh) = MDL0arr;
end



% OS scans along trajectory
i = 7;  % i=7 NOW E21 !!!
clr0 = [1 0 0]; % jet view TRUE color
clr1 = [0 0 1]; % jet view FALSE color
try
    delete(SpcftPosHnd)
    delete(SpcftVctHnd)
    delete(TargetHnd)
end
figure(1)
SpcftPosHnd = line(trj(i).trj(1,1), trj(i).trj(1,2), trj(i).trj(1,3),...
    'Marker', '.', 'MarkerSize', 40);
SpcftVctHnd = line(trj(i).trj(1:2,1), trj(i).trj(1:2,2), trj(i).trj(1:2,3),...
    'LineStyle','-','LineWidth', 3);
TargetHnd = line(trj(i).trj(1:2,1), trj(i).trj(1:2,2), trj(i).trj(1:2,3),...
    'LineStyle','-.', 'MarkerSize', 40);
Planet2ang = 0*Planet2(:,:,1);
GasMagImg = 0*Planet2(:,:,1);
GasMagImg0 = GasMagImg;
ramangARR = 3*(pi/180);%[1 3 6 10 15 20]'*(pi/180);%
phi = -20*(pi/180);%
ConeAng = 3.0*(pi/180);%3.0*(pi/180);  %Cone half angle
Dv = 100;%1;%1921;%488;%  %m/s OS delta-v HALF width
DvNum = 25;%3;%125;%
DvInt = Dv/DvNum;  %m/s
DvArr = -Dv:DvInt:Dv;  %m/s edge of Vcmp bins
DvArr(isempty(DvArr)) = 0;
DvArrCntr = (DvArr(2:end)+DvArr(1:end-1))/2;  %m/s centers of Vcmp bins
DvArrCntr(isempty(DvArrCntr)) = 0;
%TotalCnts = 0*ramangARR;
%%%%% signal figures
figure(21)
OSsigHnd = semilogy(trj(i).t, 0*PorcoSlct*(0*trj(i).t'+1));  %OS signal plot - JET RESOLVED
figure(22)
set(22,'Position',[794 474 560 420])
OSsigTOThnd = plot(trj(i).t, 0*ramangARR*trj(i).t',...
    trj(i).t, 0*trj(i).t'+1,...
    trj(i).t, 0*trj(i).t');  %OS signal plot - JET TOTAL - ALL ANGLES
axis([t0,t1,0,10*400])%1e-4,1e5])
figure(23)
set(23,'Position',[1308 205 560 420])
VcmpHnd = plot(trj(i).t, 0*ramangARR*trj(i).t');
axis([t0,t1,0,33000])
figure(76)
set(76,'Position',[1308 472 560 420])
%%%%%%%%
dspl = 1;  %full display
StripeNow = 1;
SignalMap = 0;
SignalMapCUM = 0;
SignalMapFile = 'E21H2Oosnb';%'E21H2osnb';%'E21CO2osnb';%
VcmpOnly = 0;  % only get Vcmp values, do nothing else
VcmpARR = 0*ramangARR*trj(i).t';
scan = 1;  % scan Vcomp
VcmpLoad = 1;  %load Vcmp from file - overides internal scan settings- BT 052815
VcmpFile = 'E21VcmpACTUAL.txt';%'OutputToRebecca.txt';%'E21vCmpRP1';%
VcmpPrd = 4.628;%  %s - scanning period
VcmpOffstTIME = 0.0;  %s - scanning speed time offset
AimPntTime = 0;  %s - Time used to set aimpoint below spacecraft (so desired and true aimpoints agree)
AimPntTimeInd = round(interp1q(trj(i).t, trj(i).inds', AimPntTime));
AimPoint = cnst.Re*trj(i).trj(AimPntTimeInd,:)/trj(i).rad(AimPntTimeInd,:);  % center position to "aim" scan
% VcmpScanAmp = 1000;  %m/s
% Vcmp1 = 0.5*VcmpScanAmp;  %m/s - min scanning
% Vcmp0 = -0.5*VcmpScanAmp;  %m/s - max scanning
VcmpScanAmp0 = 1000;  %m/s
%%% LOAD Vcmp File %%%
if VcmpLoad
%     VcmpLoadDat0 = load(VcmpFile);
%     VcmpLoadDat0 = VcmpLoadDat0.dat;
%     VcmpLoadDat0time = (VcmpLoadDat0(:,1)+VcmpLoadDat0(:,2))/2;
%     VcmpLoadDat0time = [trj(i).t(1); VcmpLoadDat0time; trj(i).t(end)];
%     VcmpLoadDat0 = 1000*VcmpLoadDat0(:,3);
%     VcmpLoadDat0 = [VcmpLoadDat0(1); VcmpLoadDat0; VcmpLoadDat0(end)];
%     VcmpLoadDat = interp1q(VcmpLoadDat0time, VcmpLoadDat0, trj(i).t);
%     VcmpLoadDat(isnan(VcmpLoadDat)) = 0;
%%%
    load(VcmpFile)
        load QLindH2
        inds0 = interp1q(times, (1:size(QLind,1))', trj(i).t);%min(max(trj(i).t,times(1)),times(end)));
        inds0 = round(inds0);
        QLind = interp1q((1:size(QLind,1))', QLind, inds0);
     VcmpLoadDat0 = E21VcmpACTUAL;%OutputToRebecca;%VcmpLoadDat0.dat;
     VcmpLoadDat0time = VcmpLoadDat0(:,1);%(VcmpLoadDat0(:,1)+VcmpLoadDat0(:,2))/2;
     VcmpLoadDat0time = [trj(i).t(1); VcmpLoadDat0time; trj(i).t(end)];
     VcmpLoadDat0 = VcmpLoadDat0(:,2);%1000*VcmpLoadDat0(:,3);%
     VcmpLoadDat0 = [VcmpLoadDat0(1); VcmpLoadDat0; VcmpLoadDat0(end)];
     %VcmpLoadDat = interp1q(VcmpLoadDat0time, VcmpLoadDat0, trj(i).t);
        %%% avoids interpolating right between Vcmp scans, can give false signal peaks %%%
        indVcmp = interp1q(VcmpLoadDat0time, (1:numel(VcmpLoadDat0))', trj(i).t);
        cndVcmp = ((VcmpLoadDat0(ceil(indVcmp))-VcmpLoadDat0(floor(indVcmp)))<0);
        indVcmp(cndVcmp) = floor(indVcmp(cndVcmp));
        VcmpLoadDat = interp1q((1:numel(VcmpLoadDat0))', VcmpLoadDat0, indVcmp);
        %%%%%%%%%%%%%%%%%%%%%
     VcmpLoadDat(isnan(VcmpLoadDat)) = 0;
end
%%%%%%%%%%%%%%%%%%%%%%

% load SignalMapDat  
%%%time shift something messed up - fix later - BT 011716 !!!!!
load(SignalMapFile)
SignalMapDat0time = SignalMapDat(:,1);
SignalMapDat0 = SignalMapDat(:,2);
SignalMapDat0 = [SignalMapDat0(1); SignalMapDat0; SignalMapDat0(end)];
SignalMapDat = interp1q(SignalMapDat0time, SignalMapDat0, trj(i).t);
dtmin = 0.15;%min(SignalMapDat0time(2:end)-SignalMapDat0time(1:end-1));
ind = floor(interp1q(SignalMapDat0time, (1:numel(SignalMapDat0time))', trj(i).t));
cnd = ((trj(i).t - SignalMapDat0time(ind)) < 2*dtmin);
SignalMapDat(cnd==0) = 0;
%plot(trj(i).t,SignalMapDat,trj(i).t,10*cnd)
%set(OSsigTOThnd(end), 'XData', SignalMapDat(:,1), 'YData', SignalMapDat(:,2))
set(OSsigTOThnd(end), 'XData', trj(i).t-0.7*dtmin, 'YData', SignalMapDat)

%%%% specify feature times & indices %%%%
%peakTimes = [-1.86, -1.11, -0.66, -0.22, 0.59, 1.68, 2.2, 2.49]; %H2
%peakTimes = [-1.86, -0.22, 0.59, 1.68]; %H2
peakTimes = [-12.4, -8.02, 0.76, 2.12, 8.382, 19.41]; %CO2
iiiarr2 = iiiarr;
%iiiarr2 = round(interp1q(trj(i).t,iiiarr2',peakTimes'))';
%%%%%%%



%% MAIN LOOPS %%%
%plot(trj(i).t,50*OSNBsens*JetSigOStot, '-', trj(i).t, 10*SignalMapDat, '-')
% ramangARR = 3*(pi/180);
% phi = -20*(pi/180);%

 
%  ramangARR = 7.5*(pi/180);
%  phi = -43*(pi/180);

%  ramangARR = 23*(pi/180);
%  phi = 15*(pi/180);

%  ramangARR = 7*(pi/180);
%  phi = 20*(pi/180);

 ramangARR = 3*(pi/180);
 phi = -20*(pi/180);

    %OSNBsens = 4*3e-16;%2000*3e-16;%
    
 %set(OSsigTOThnd(end), 'XData', trj(i).t, 'YData', SignalMapDat/10000)
 %axis([t0,t1,0,1000])
 GasMagImg0 = 0*GasMagImg0;
 GasMagImg = 0*GasMagImg;
tic
    OSNBsection2b%_alt%_ISOTROP
toc

end