%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: simulateTraj.m
% Author: J. Kristoffer Hultgren
%
% Description:
%
% Simulates a sweep during a trajectory
%
% Input:  -
%
% Output: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function simulateTraj()

LP_Constants;

dist = 1.3;         % Change to that of use
heli = '1.3';       % Change to that of use

if strcmp(heli,'1.3')
    %disp('Run orbit.dat in GUI for DSMC and MHD models, save model output in simu/, then press any key to continue.')
    %pause
    
    % Refine the data file
    %!sed -i '1,3d' simu/MHD_CG_1.3_au.dat
    %!sed -i '1,3d' simu/DSMC_CG_1.3_au.dat
    
    % For compatibility
    %!mv simu/MHD_CG_1.3_au.dat simu/MHD13.dat
    %!mv simu/DSMC_CG_1.3_au.dat simu/DSMC13.dat

    % Load file
    load('simu/DSMC13.dat');
    
    % Controls size of data set (occationally the data set recieved for DSMC
    %                               was missing the last line)
    nCompare = numel(DSMC13(:,1));
    half = floor(nCompare/2);
    if 2*half == nCompare
        %DSMC13 = [DSMC13;DSMC13(nCompare,:)];
    end

    % Load file
    load('simu/MHD13.dat');

    % Separate data
    ni  = MHD13(:,5).*1e-6;
    nhi = zeros(size(ni));
    vi  = sqrt(MHD13(:,6).^2 + MHD13(:,7).^2 + MHD13(:,8).^2).*10^-3;
    vhi = zeros(size(vi));
    
    data = MHD13;
     
    MT = MHD13(:,12)./((MHD13(:,5) + DSMC13(:,9)).*CO.Kb1).*CO.Kb2;
    nT = DSMC13(:,8).*CO.Kb2;
    weightT = (MHD13(:,5).*MT+DSMC13(:,9).*nT)./(MHD13(:,5)+DSMC13(:,9));
    
    Ti = nT;
    Te = Ti;
    
    Bx = MHD13(:,9).*1e9;
    By = MHD13(:,10).*1e9;
    Bz = MHD13(:,11).*1e9;
    B  = (Bx.^2 + By.^2 + Bz.^2).^(0.5);
        
end

if strcmp(heli,'2.0')
    
    %disp('Run orbit.dat in GUI for DSMC and MHD models, save model output in simu/, then press any key to continue.')
    %pause
    
    % Refine the data file
    %!sed -i '1,3d' simu/MHD_CG_2.0_au.dat
    %!sed -i '1,3d' simu/DSMC_CG_2.0_au.dat
    
    % For compatibility
    %!mv simu/MHD_CG_2.0_au.dat simu/MHD20.dat
    %!mv simu/DSMC_CG_2.0_au.dat simu/DSMC20.dat

    % Load file
    load('simu/DSMC20.dat');
    
    % Controls size of data set (occationally the data set recieved for DSMC
    %                               was missing the last line)
    nCompare = numel(DSMC20(:,1));
    half = floor(nCompare/2);
    if 2*half == nCompare
        %DSMC20 = [DSMC20;DSMC20(nCompare,:)];
    end

    % Load file
    load('simu/MHD20.dat');

    % Separate data
    ni  = MHD20(:,5).*1e-6;
    nhi = zeros(size(ni));
    vi  = sqrt(MHD20(:,6).^2 + MHD20(:,7).^2 + MHD20(:,8).^2).*10^-3;
    vhi = zeros(size(vi));
    
    data = MHD20;
     
    MT = MHD20(:,12)./((MHD20(:,5) + MHD20(:,9)).*CO.Kb1).*CO.Kb2;
    nT = DSMC20(:,8).*CO.Kb2;
    weightT = (MHD20(:,5).*MT+DSMC20(:,9).*nT)./(MHD20(:,5)+DSMC20(:,9));
    
    Ti = nT;
    Te = Ti;
    
    Bx = MHD20(:,9).*1e9;
    By = MHD20(:,10).*1e9;
    Bz = MHD20(:,11).*1e9;
    B  = (Bx.^2 + By.^2 + Bz.^2).^(0.5);
    
end

if strcmp(heli,'2.7')
    
    %disp('Running Hybrid model...')
    
    % Run model
    %!./ICES.pl orbit.dat CG_2.7_au simu HYB
        
    %disp('Run orbit.dat in GUI for DSMC model, save model output in simu/, then press any key to continue.')
    %pause
    
    % Refines the data files
    %!sed -i '1,34d' simu/Hybrid_CG_2.7_au.dat
    %!sed -i '1,3d' simu/DSMC_CG_2.7_au.dat
    
    % For compatibility
    %!mv simu/Hybrid_CG_2.7_au.dat simu/Hybrid27.dat
    %!mv simu/DSMC_CG_2.7_au.dat simu/DSMC27.dat

    % Load file and add additional rows for compatibility
    load('simu/DSMC27.dat');
    
    % Controls size of data set (occationally the data set recieved for DSMC
    %                               was missing the last line)
    nCompare = numel(DSMC27(:,1));
    half = floor(nCompare/2);
    if 2*half == nCompare
     %   DSMC27 = [DSMC27;DSMC27(nCompare,:)];
    end

    % Load file
    load('simu/Hybrid27.dat');

    % Separate data
    ni  = Hybrid27(:,11).*1e-6;
    nhi = Hybrid27(:,15).*1e-6;
    vi  = sqrt(Hybrid27(:,12).^2 + Hybrid27(:,13).^2 + Hybrid27(:,14).^2).*10^-3;
    vhi = sqrt(Hybrid27(:,16).^2 + Hybrid27(:,17).^2 + Hybrid27(:,18).^2).*10^-3;
    Tn  = DSMC27(:,8);
    data = DSMC27;
    
    Te = zeros(size(Tn));
    Ti = Te;
    Te(:) = 1;     % Assumed value [eV]
    Ti(:) = 0.1;   % Assumed value [eV]
    
    Bx = Hybrid27(:,5).*1e9;
    By = Hybrid27(:,6).*1e9;
    Bz = Hybrid27(:,7).*1e9;
    B  = (Bx.^2 + By.^2 + Bz.^2).^(0.5);
    
end

if strcmp(heli,'3.3')
    
    %disp('Running Hybrid model...')
    
    % Run model
    %!./ICES.pl orbit.dat CG_3.3_au simu HYB
    
    %disp('Run orbit.dat in GUI for DSMC model, save model output in simu/, then press any key to continue.')
    %pause

    % Refines the data files
    %!sed -i '1,34d' simu/Hybrid_CG_3.3_au.dat
    %!sed -i '1,3d' simu/DSMC_CG_3.3_au.dat
    
    % For compatibility
    %!mv simu/Hybrid_CG_3.3_au.dat simu/Hybrid33.dat
    %!mv simu/DSMC_CG_3.3_au.dat simu/DSMC33.dat

    % Load file and add additional rows for compatibility
    load('simu/DSMC33.dat');
    
    % Controls size of data set (occationally the data set recieved for DSMC
    %                               was missing the last line)
    nCompare = numel(DSMC33(:,1));
    half = floor(nCompare/2);
    if 2*half == nCompare
        %DSMC33 = [DSMC33;DSMC33(nCompare,:)];
    end

    % Load file
    load('simu/Hybrid33.dat');

    % Separate data
    ni  = Hybrid33(:,11).*1e-6;
    nhi = Hybrid33(:,15).*1e-6;
    vi  = sqrt(Hybrid33(:,12).^2 + Hybrid33(:,13).^2 + Hybrid33(:,14).^2).*10^-3;
    vhi = sqrt(Hybrid33(:,16).^2 + Hybrid33(:,17).^2 + Hybrid33(:,18).^2).*10^-3;
    Tn  = DSMC33(:,8);
    data = DSMC33;
    
    Te = zeros(size(Tn));
    Ti = Te;
    Te(:) = 1;     % Assumed value [eV]
    Ti(:) = 0.1;   % Assumed value [eV]
    
    Bx = Hybrid33(:,5).*1e9;
    By = Hybrid33(:,6).*1e9;
    Bz = Hybrid33(:,7).*1e9;
    B  = (Bx.^2 + By.^2 + Bz.^2).^(0.5);
    
end

Vb = 0.25;           % Bias
nVb = max(size(Vb));        % # bias steps

disp('Done!')
disp('Running instrument model...')

% Calculate currents
for i = 1:max(size(ni))
    for k = 1:nVb
        temp = inst(dist,ni(i),nhi(i),Te(i),Ti(i),vi(i),vhi(i),Vb(k));
        currout(i,k) = temp(1);
    end
    disp('Percent done:')
    disp(i*100/max(size(ni)))
end

if nVb == 1
    % Calculates currents
    for i = 1:max(size(ni))
        currout2(i,:) = inst(dist,ni(i,1),nhi(i,1),Te(i,1),Ti(i,1),vi(i,1),vhi(i,1),Vb);
    end
end

disp('Done!')
disp('Rearranging and plotting data...')

%Lower resolution to that of the measurement ranges (simulate ADC):
span    = 10e-6;                            %Gain spanning (A)
deltaI  = span/2^16;                        %Gain spanning / #bits 
currout = deltaI.*round(currout./deltaI);
if nVb == 1
    currout2 = deltaI.*round(currout2./deltaI);
    currout2 = currout2.*1e9;
end

l = linspace(data(1,1),data(max(size(data)),1),max(size(data))).*1e-3;
currout = currout.*1e9;

%Plot sweep during time of flight
%titel = 'I_{p} (nA)';
%maxx = max(max(currout));
%meann = mean(mean(currout));
%if (maxx - meann) > (3*maxx/4)
%    currout = log10(abs(currout));
%    titel = 'I_{p} log_{10}(nA)';
%end

%figure('name','Timeplot Sweep')
%pcolor(l,Vb,currout')
%shading interp
%colorbar
%title(titel)
%xlabel('Distance (km)')
%ylabel('V_{b} (V)')

%Plot parameters used
figure('name','Parameters')
subplot(3,1,1), plot(l,Bx,'b')
title('Plasma Parameters')
hold on
subplot(3,1,1), plot(l,By,'g')
subplot(3,1,1), plot(l,Bz,'r')
ylabel('B-field (nT)')
legend('B_{x}','B_{y}','B_{z}')
subplot(3,1,2), plot(l,Ti)
ylabel('Fluid temperature (eV)')
subplot(3,1,3), plot(l,log10(ni+nhi))
hold on
%ylabel('Fluid density log_{10}(cm^{-3})')
plot(l,log10(currout2(:,1)),'g')
%plot(l,log10(currout2(:,2)),'r')
%plot(l,log10(currout2(:,3)),'g')
legend('Fluid density log_{10}(cm^{-3})','I_{p} log_{10}(nA)')
%ylabel('I_{p} log_{10}(nA)')
xlabel('Distance (km)')


disp('Done!')
disp('Finished!')