% load test2 FitPrmsOpt JetLines
%fct = 1;%dat(3).CalFct;  %for test 1 and 2 only !!!
%figure(1)

%MDLPrv = 0*MDL;draw
SingleJet = 1;


SetExact = 1;  % true setting screws up optimization, fix later
%SetExact = 0;  % true setting screws up optimization, fix later
if SetExact
    if LoadPorco
        if ~exist('PorcoDat', 'var')
            load PorcoDat % Lat, Wlon, TiltZen, TiltAz, NormSight
            PorcoDat0 = PorcoDat;
        else
            PorcoDat = PorcoDat0;
        end
    PorcoDat(:,1:4) = (pi/180)*PorcoDat0(:,1:4);
    crds0 = 0*[PorcoDat(:,1:3),PorcoDat(:,1:3)];
    lon0 = 0*(pi/180);
    az0 = 180*(pi/180);
    PorcoDat(:,2) = -PorcoDat(:,2) + lon0;
    PorcoDat(:,4) = PorcoDat(:,4) + az0;

    tilt = PorcoDat(:,3);
    azimuth = PorcoDat(:,4);
    lat_polar = (pi/2) - PorcoDat(:,1);
    longitude = PorcoDat(:,2);
    % jet on top
    crds0(:,1) = 0; %x
    crds0(:,2) = 0; %y
    crds0(:,3) = 1; %z
    crds0(:,4) = 0;
    crds0(:,5) = 0;
    crds0(:,6) = 1;
    % tilt jet (about y-axis)
    vx = crds0(:,4).*cos(tilt) - crds0(:,6).*sin(tilt);
    vy = crds0(:,5);
    vz = crds0(:,4).*sin(tilt) + crds0(:,6).*cos(tilt);
    crds0(:,4:6) = [vx vy vz];
    % rotate jet to specified azimuth (i.e. about z-axis)
    vx = crds0(:,4).*cos(azimuth) + crds0(:,5).*sin(azimuth);
    vy = -crds0(:,4).*sin(azimuth) + crds0(:,5).*cos(azimuth);
    vz = crds0(:,6);
    crds0(:,4:6) = [vx vy vz];
    %% rotate source to specified latitude (about y-axis)
    x = crds0(:,1).*cos(lat_polar) - crds0(:,3).*sin(lat_polar);
    y = crds0(:,2);
    z = crds0(:,1).*sin(lat_polar) + crds0(:,3).*cos(lat_polar);
    vx = crds0(:,4).*cos(lat_polar) - crds0(:,6).*sin(lat_polar);
    vy = crds0(:,5);
    vz = crds0(:,4).*sin(lat_polar) + crds0(:,6).*cos(lat_polar);
    crds0(:,1:6) = [x y z vx vy vz];
    %% rotate source to specified longitude (i.e. about z-axis)
    x = crds0(:,1).*cos(-longitude) + crds0(:,2).*sin(-longitude);
    y = -crds0(:,1).*sin(-longitude) + crds0(:,2).*cos(-longitude);
    z = crds0(:,3);
    vx = crds0(:,4).*cos(-longitude) + crds0(:,5).*sin(-longitude);
    vy = -crds0(:,4).*sin(-longitude) + crds0(:,5).*cos(-longitude);
    vz = crds0(:,6);
    crds0(:,1:6) = [x y z vx vy vz];
    crds0(:,1:3) = cnst.Re*crds0(:,1:3);
    %end
    
    FitPrms = zeros(size(crds0,1),14+3);
    FitPrms(:,11:13) = crds0(:,4:6);  %src jet direction (normalized)
    FitPrms(:,8) = Mach;  % Mach number
    %FitPrms(:,5) = ii;  % which TS
    if ~EnhanceON
        FitPrms(:,10) = SrcRt0;%3*PorcoDat(:,5)*SrcRt0;%   %SRC RATE  PorcoDat(:,5)*SrcRt0;
    else
        FitPrms(:,10) = AvgEnh*SrcRt0;   %SRC RATE
    end
    FitPrms(:,14:16) = crds0(:,1:3);  %src position
    
    else
    
    FitPrms0 = zeros(0,14+3);
    indVl0 = 0;
    TgrT = [];
    for ii = 1:length(TgrStrps)
        %ii = 1;
        crds0 = TgrStrps(ii).crds;
        crds0NrmMG = sqrt(dot(crds0,crds0,2));
        crds0Nrm = crds0 ./ [crds0NrmMG, crds0NrmMG, crds0NrmMG];
        crds0 = [crds0, crds0Nrm];
        %%%%
        TgrT = [TgrT; TgrStrps(ii).T];
        
        tilt = 0;%PorcoDat(:,3);
        azimuth = 0;%PorcoDat(:,4);
        % tilt jet (about y-axis)
        vx = crds0(:,4).*cos(tilt) - crds0(:,6).*sin(tilt);
        vy = crds0(:,5);
        vz = crds0(:,4).*sin(tilt) + crds0(:,6).*cos(tilt);
        crds0(:,4:6) = [vx vy vz];
        % rotate jet to specified azimuth (i.e. about z-axis)
        vx = crds0(:,4).*cos(azimuth) + crds0(:,5).*sin(azimuth);
        vy = -crds0(:,4).*sin(azimuth) + crds0(:,5).*cos(azimuth);
        vz = crds0(:,6);
        crds0(:,4:6) = [vx vy vz];
        
        
        %     crds0(:,1:3) = cnst.Re*crds0(:,1:3);
        %     %crds0(:,4:6) = 0.25*JetLen*crds0(:,4:6);
        FitPrms = zeros(size(crds0,1),14+3);
        FitPrms(:,11:13) = crds0(:,4:6);
        FitPrms(:,8) = Mach;  % Mach number
        FitPrms(:,5) = ii;  % which TS
        %FitPrms(:,10) = SrcRt0nrm*TgrStrps(ii).src;
        if ~EnhanceON
            FitPrms(:,10) = SrcRt0nrm*TgrStrps(ii).src;%   %SRC RATE  PorcoDat(:,5)*SrcRt0;
        else
            indVls = indVl0 + (1:numel(TgrStrps(ii).src));
            FitPrms(:,10) = SrcRt0nrm*AvgEnh(indVls).*TgrStrps(ii).src;   %SRC RATE
        end
        FitPrms(:,14:16) = crds0(:,1:3);  %src position
        
        FitPrms0 = [FitPrms0; FitPrms];
        indVl0 = indVl0 + numel(TgrStrps(ii).src);
    end
    FitPrms = FitPrms0;
    
    
    PorcoSlct = 0*(1:size(FitPrms,1))' + 1;  % select all
    
    end
    

    
%     FitPrms0 = [FitPrms0; FitPrms];
%     FitPrms = FitPrms0;

%%%% SOURCE COLORS %%%%%
if exist('AvgEnh')
    if std(AvgEnh)
        EnhanceON = 1;
    end
end
PorcoDsg = 0*FitPrms(:,10);  % designates tiger stripe source
if 1%TgrStrpColor
    %if LoadPorco% && EnhanceON
        PorcoDsg([1:24,94:95]) = 1;  %Cairo
        PorcoDsg([25:59,96:97]) = 2; %Baghdad
        PorcoDsg([60:88]) = 3;       %Damascas
        PorcoDsg([89:93,98]) = 4;    %Alexandria
        PorcoClr = 0*[PorcoDsg,PorcoDsg,PorcoDsg];
        PorcoClr(PorcoDsg==1,:) = (0*PorcoClr(PorcoDsg==1,1)+1)*[1 0 0];
        PorcoClr(PorcoDsg==2,:) = (0*PorcoClr(PorcoDsg==2,1)+1)*[0 0 1];
        PorcoClr(PorcoDsg==3,:) = (0*PorcoClr(PorcoDsg==3,1)+1)*[0 1 0];
        PorcoClr(PorcoDsg==4,:) = (0*PorcoClr(PorcoDsg==4,1)+1)*[1 0 1];
        
        %fvl = 0.8*AvgEnh/6; % your float
        if EnhanceON
            colormap('parula')
            cm = colormap; % returns the current color map
            
            if LoadPorco%EnhanceON
                fvl = 0.8*AvgEnh/6; % your float
            else
                fvl = (TgrT - min(TgrT(:)))/(max(TgrT(:)) - min(TgrT(:))); % your float
            end
            colorID = sum(fvl*(0*[0:length(cm(:,1))] + 1) >...
                (0*(1:size(FitPrms,1))+1)'*[0:1/length(cm(:,1)):1],2);
            colorID = max(1, colorID);
            colorID = min(64, colorID);
            PorcoClr = cm(colorID, :); % returns your color
        end
    %end
    
end
%%%%%%%%%%



%% Mach numbers & Source rates

m = 18*1.67e-27;  %kh molecule mass
k = 1.38e-23;  %J/K  Boltzmann constant
vth = sqrt(8*k*Tsrc/pi/m);
prfcts = FitPrms(:,10);
Mvls = FitPrms(:,8);
RamFct = (((2/sqrt(pi)) - (1./Mvls)).*exp(-Mvls.^2) +...
    (2*Mvls + (1./Mvls)).*(1+erf(Mvls)));
RamFct(isnan(RamFct)) = 2*(2/sqrt(pi));  % M=0 condition
%F0degThrstrTST = exp(-S0tst.^2) +  sqrt(pi)*S0tst.*(1 + erf(S0tst));
    %RamFct = RamFct.*Mvls;
SrcRts = vth*prfcts.*RamFct;  %molecules/sec Dong et al
% Mvls(SrcRts==0) = [];
% FitPrms(SrcRts==0,:) = [];
% %JetLines(SrcRts==0) = [];
% SrcRts(SrcRts==0) = [];
% %[Mvls, SrcRts*m]
TotalSource = sum(SrcRts)*m



end


%% go through all jets
%if 0
MDL0arr = zeros(size(FitPrms,1),numel(trj(1).t),flbylen);
for rr = 1:size(FitPrms,1)

    if PorcoSlct(rr)
%[cnt, Dt, DevVl, swvl, ii, jj, kk, M, thtJet, prfct, aimVctNrm, qqq];
             
% unpackage
%cnt = FitPrms(rr,1);
%Dt = FitPrms(rr,2);
%DevVl = FitPrms(rr,3);
%swvl = FitPrms(rr,4);
ii = FitPrms(rr,5);  % which tiger stripe
jj = FitPrms(rr,6);  % where on tiger stripe
%kk = FitPrms(rr,7);
M = FitPrms(rr,8);
%thtJet = FitPrms(rr,9);
prfct = FitPrms(rr,10);
aimVctNrm = FitPrms(rr,11:13);
%qqq = FitPrms(rr,14);

% set exact source positions
if ~SetExact
    crds = TgrStrps(ii).crds(jj,:);
else
    %aimVctNrm = FitPrms(rr,11:13);%aimVctNrmVls(rr,:);
    crds = FitPrms(rr,14:16);%crds0(rr,:);
end


% show result
PlumeModelCalc2

% scale by TrueAnomFct
MDL = ((0*(1:size(MDL,1))+1)'*TrueAnomFct).*MDL;

% add previous 'background' fit
MDL0(rr).MDL = MDL;  % record individual profile
MDL0arr(rr,:,:) = MDL;
MDL = MDL + MDLPrv;
MDLPrv = MDL;  % Record model output
            
%rr
end
%end


%% constructs jet lines and update display
if EnhanceON
    JetLen = 2*AvgEnh(rr)*0.05e5;  %m Jet display length
else
    JetLen = 1*0.25*1.0e5;  %m Jet display length
end
if DrawJets
if LoadPorco
    % all jets
    try %~exist('JetHnd')
        set(JetHnd(rr),...
            'XData', [crds(1) crds(1)+JetLen*aimVctNrm(1)],...
            'YData', [crds(2) crds(2)+JetLen*aimVctNrm(2)],...
            'ZData', [crds(3) crds(3)+JetLen*aimVctNrm(3)],...
            'Parent',AxHnd,...
            'LineWidth',3, 'Color', PorcoClr(rr,:),...
            'Visible', 'on');
    catch
        JetHnd(rr) = line('XData', [crds(1) crds(1)+JetLen*aimVctNrm(1)],...
            'YData', [crds(2) crds(2)+JetLen*aimVctNrm(2)],...
            'ZData', [crds(3) crds(3)+JetLen*aimVctNrm(3)],...
            'Parent',AxHnd,...
            'LineWidth',3, 'Color', PorcoClr(rr,:),...
            'Visible', 'on');
    end


% % enhanced jets
% JetLen2 = 1*JetLen;
% JetHnd2(rr) = line('XData', [crds(1) crds(1)+JetLen2*aimVctNrm(1)],...
%     'YData', [crds(2) crds(2)+JetLen2*aimVctNrm(2)],...
%     'ZData', [crds(3) crds(3)+JetLen2*aimVctNrm(3)],...
%     'Parent',AxHnd,...
%     'LineWidth',1*2, 'Color', PorcoClr(rr,:),...
%     'Visible', 'off');


end


% source dots
try
    set(SrcDotHnds(rr),'Color',PorcoClr(rr,:))
catch
    SrcDotHnds(rr) = line('XData', [crds(1) + 0.01*cnst.Re*aimVctNrm(1)],...
        'YData', [crds(2) + 0.01*cnst.Re*aimVctNrm(2)],...
        'ZData', [crds(3) + 0.01*cnst.Re*aimVctNrm(3)],...
        'Parent',AxHnd,...
        'Color', PorcoClr(rr,:),...
        'Marker', '.',...
        'MarkerSize', 1.4*40);
end


end

end
%end
    

            


%%  update MDL display (optionaly)
if Plt
for j = AnalysisTrj
    
    %Juergen Files
    if 0
        if j == 3 %E7 data
            fnm = 'CDA_hrd_2009_306_E7_1.7_micron.ascii';
            %fnm = 'CDA_hrd_2009_306_E7_3.0_micron.ascii';
            tcol = 2;
            DustScl = 40e2;
        elseif j == 5 %E14 model
            fnm = 'E14Profiles_Model42_ascii.out';
            tcol = 1;
            DustScl = 5e2;
        elseif j == 4 %E17 model
            fnm = 'E17Profiles_Model42_ascii.out';
            tcol = 1;
            DustScl = 5e2;
        else  %DEFAULT - No Juergen file available
            fnm = 'E14Profiles_Model42_ascii.out';
            tcol = 1;
            DustScl = 5e2;
        end
        DustDat = importdata(fnm,' ');
        DustTime = DustDat.data(:,tcol);
        DustDns = DustDat.data(:,3);
        DustDns = interp1q(DustTime, DustDns, dat(j).t);
        amp = DustScl*DustDns;
    end
    
    
    
    % set fit surface Z
    % cancel instrument CalFct and apply by display scale factor
    amp = (dat(j).sclfct/dat(j).CalFct)*MDL(:,j); % MDL is CD for UVIS, density for INMS !!!
    %amp = (dat(j).sclfct/1)*MDL(:,j);%(5e-17)*MDL(:,j);%
    sigsurf_MDL(:,:,3,j) = sigsurf_MDL0(:,:,3,j) - ((N1arr-1)/N1)*amp';
    
    % display fit line
    if DrawFitLine
        clr = [0 0 0]+0;
        line('XData', sigsurf_MDL(1,:,1,j),...
            'YData', sigsurf_MDL(1,:,2,j),...
            'ZData', sigsurf_MDL(1,:,3,j)-amp',...
            'LineWidth', 4, 'Color',clr,...
            'Parent',AxHnd)
    end
    
    
%         line('XData', trj(1).trj(:,1),...
%             'YData', trj(1).trj(:,2),...
%             'ZData', trj(1).trj(:,3)-amp,...
%             'LineWidth', 4, 'Color',clr)%,...
%             %'Parent',AxHnd)

end
end
% option to show jet vector as it's being adjusted
% if dspCurJet
%     set(JetHnd(end),'Visible','on')
% else
%     set(JetHnd(end),'Visible','off')
% end

%%
if Plt
    PlottingInstructions2
    pause(0.3)
end