%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: simulate.m
% Author: J. Kristoffer Hultgren
%
% Description:
%
% Simulates LAP measurements and outputs the resulting current/voltage.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function simulate()

%Asking for orbital elements
%dist = input('What heliocentric distance do you want to model? (1.3/2.0/2.7/3.3): ');
dist = 1.3;
heli='1.3';% = num2str(dist);
r = 5;% = input('What is the radius of the orbit (km)? ');
r = r*1000;
inc=pi/2;% = input('What is the inclination of the orbit (rad)? ');
omega=3*pi/2;% = input('What is the right ascension (rad)? ');

disp('Initializing...')

LP_Constants;        %Defines physical constants

elem = draworbit(r,inc,omega);  %Draws the orbit, output the x,y,z-coordinates as a text file, and returns the orbital elements as a matrix

v = sqrt(CO.G*CO.M./r);   %Speed in orbit (m/s)

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
    
    % Refines the data files
    %!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
    
    % Plot data in orbit
    plotDSMCPolarOrbit(DSMC13,elem,v);
    
    % Load file
    load('simu/MHD13.dat');
    
    % Plot data in orbit
    plotMHDPolarOrbit(MHD13,elem); %#ok<NODEF>
    plotMHDTime(MHD13,r);
    
    % Separate data
    ni  = MHD13(:,5).*1e-6;
    nhi = zeros(size(ni));
    vi  = sqrt(MHD13(:,6).^2 + MHD13(:,7).^2 + MHD13(:,8).^2).*1e-3;
    vhi = zeros(size(vi));
    
    data = MHD13;
    
    nT  = DSMC13(:,8).*CO.Kb2;
    MT = MHD13(:,12)./((MHD13(:,5) + DSMC13(:,9)).*CO.Kb1).*CO.Kb2;
    weightT = (MHD13(:,5).*MT+DSMC13(:,9).*nT)./(MHD13(:,5)+DSMC13(:,9));
    
    Ti = nT;
    Te = Ti;
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
    
    % Refines 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
        
    % Plot data in orbit
    plotDSMCPolarOrbit(DSMC20,elem,v);
    
    % Load file
    load('simu/MHD20.dat');
    
    % Plot data in orbit
    plotMHDPolarOrbit(MHD20,elem); %#ok<NODEF>
    plotMHDTime(MHD20,r);
    
    % Separate data
    ni  = MHD20(:,5).*1e-6;
    nhi = zeros(size(ni));
    vi  = sqrt(MHD20(:,6).^2 + MHD20(:,7).^2 + MHD20(:,8).^2).*1e-3;
    vhi = zeros(size(vi));
  
    data = MHD20;
    
    nT = DSMC20(:,8).*CO.Kb2;
    MT = MHD20(:,12)./((MHD20(:,5) + DSMC20(:,9)).*CO.Kb1).*CO.Kb2;
    weightT = (MHD20(:,5).*MT+DSMC20(:,9).*nT)./(MHD20(:,5)+DSMC20(:,9));
    
    Ti = nT;
    Te = Ti;
        
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
    
    % Plot data in orbit
    plotDSMCPolarOrbit(DSMC27,elem,v);
    
    % Load file
    load('simu/Hybrid27.dat');
    
    % Plot data in orbit
    plotHybridPolarOrbit(Hybrid27,elem,v); %#ok<NODEF>
    plotHybridTime(Hybrid27,r);
    
    % 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 = Hybrid27;
  
    Te = zeros(size(Tn));
    Ti = Te;
    Te(:) = 1;     % Assumed value [eV]
    Ti(:) = 0.1;   % Assumed value [eV]
    
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
    
    % Plot data in orbit
    plotDSMCPolarOrbit(DSMC33,elem,v);
    
    % Load file
    load('simu/Hybrid33.dat');
    
    % Plot data in orbit
    plotHybridPolarOrbit(Hybrid33,elem,v); %#ok<NODEF>
    plotHybridTime(Hybrid33,r);
    
    % Separate data
    ni  = Hybrid33(:,11).*1e-6;
    nhi = Hybrid33(:,15).*1e-6;
    vi  = sqrt(Hybrid33(:,12).^2 + Hybrid33(:,13).^2 + Hybrid33(:,14).^2).*1e-3;
    vhi = sqrt(Hybrid33(:,16).^2 + Hybrid33(:,17).^2 + Hybrid33(:,18).^2).*1e-3;
    Tn  = DSMC33(:,8);
    data = Hybrid33;
    
    Te = zeros(size(Tn));
    Ti = Te;
    Te(:) = 1;     % Assumed value [eV]
    Ti(:) = 0.1;   % Assumed value [eV]
    
end

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

Vb = 1;     % Bias voltage

currout = zeros(max(size(ni)),3);

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

disp('Done!')
disp('Rearranges and plots 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);

currout = currout.*1e9;

ne = ni + nhi;      % Electron density

plotCurrents(data,currout(:,1),currout(:,2),currout(:,3),r,ne,Te);

disp('Done!')

disp('Speed in orbit (m/s):')
disp(v)

disp('Finished!')




